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BACKGROUND OF THE INVENTION 

a . Field of Invention 

The invention relates generally to computerized 
15 techniques for processing data obtained from radar to track 
multiple discrete objects. 

b. Description of the Background 

There are many situations where the courses of multiple 
objects in a region must be tracked. Typically, radar is used 

2 0 to scan the region and generate discrete images or "snapshots" 
based on sets of returns or observations. In some types of 
tracking systems, all the returns from any one object are 
represented in an image as a single point unrelated to the 
shape or size of the objects. "Tracking" is the process of 

25 identifying a sequence of points from a respective seq^>e^lce of 
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the images that represents the motion of an object. The 
tracking problem is difficult when there are multiple closely 
spaced objects because the objects can change speed and 
direction rapidly and move into and out of the line of sight 
5 for other objects. The problem is exacerbated because each 
set of returns may result from noise as well as echoes from 
the actual objects. The returns resulting from the noise are 
also called false positives. Likewise, the radar will not 
detect all echoes from the actual objects and this phenomena 

10 is called a false negative or "missed detect" error. For 
tracking airborne objects, a large distance between the radar 
and the objects diminishes the signal to noise ratio so the 
number of false positives and false negatives can be high. 
For robotic applications, the power of the radar is low and as 

15 a result, the signal to noise ratio can also be low and the 
number of false positives and false negatives high. 

In view of the proximity of the objects to one another, 
varied motion of the objects and false positives and false 
negatives, multiple sequential images should be analyzed 

20 collectively to obtain enough information to properly assign 
the points to the proper tracks. Naturally, the larger the 
number of images that are analyzed, the greater the amount of 
information that must be processed. 
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While identifying the track of an object, a kinematic 
model may be generated describing the location, velocity and 
acceleration of the object. Such a model provides the means 
by which the future motion of the object can be predicted. 
5 Based upon such a prediction, appropriate action may be 
initiated. For example, in a military application there is a 
need to track multiple enemy aircraft or missiles in a region 
to predict their objective, plan responses and intercept them. 
Alternatively, in a commercial air traffic control application 

10 there is a need to track multiple commercial aircraft around 
an airport to predict their future courses and avoid 
collision. Further, these and other applications, such as 
robotic applications, may use radar, sonar, infrared or other 
object detecting radiation bandwidths for tracking objects. 

15 In particular, in robotic applications reflected radiation can 
be used to track a single object which moves relative to the 
robot (or vice versa) so the robot can work on the object. 

Consider the very simple example of two objects being 
tracked and no false positives or false negatives. The radar, 

20 after scanning at time t^^, reports objects at two locations in 
a first observation set. That is, the radar returns a set of 
two observations {o^^, } - At time the radar returns a 
similar set of two observations {021, O22) from a second 
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observation set. Suppose from prior processing that track 
data for two tracks and T2 includes the locations at to of 
two objects. Track may be extended through the points in 
the two sets of observations in any of four ways, as may track 
5 T2. The possible extensions of can be described as: {T^, 

^11 / ^22 }/ {^1/ O12/ O21} and {Tj, Oj_2, O22} - Tracks 
can likewise be extended from T2 in four possible ways, 
including {T2, Oj^2f O22} - Figure 1 illustrates these five (out 
of eight) possible tracks (to simplify the problem for 

10 purposes of explanation) , The five track extensions are 
labeled h^^, h:^2i -^23/ ^i^/ and wherein h;^^ is derived from [T^, 
Oil, is derived from {T^, o^i, O22}/ is derived from 

{^1/ 0^2/ O21) , hi4 is derived from {T^, 0^2/ O22}/ and ^22 is 
derived from {T2, 0^, O21] - The problem of determining which 

15 such track extensions are the most likely or optimal is 
hereinafter known as the assignment problem. 

It is known from prior art to determine a figure of merit 
or cost for assigning each of the points in the images to a 
track. The figure of merit or cost is based on the likelihood 

20 that the point is actually part of the track. For example, 
the figure of merit or cost may be based on the distance from 
the point to an extrapolation of the track. Figure 1 
illustrates costs c52i and 622 fo^ hypothetical extension and 
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modeled target characteristics. The function to calculate the 
cost will normally incorporate detailed characteristics of the 
sensor, such as probability of measurement error, and track 
characteristics, such as likelihood of track maneuver. 
5 Figure 2 illustrates a two by two by two matrix, c, that 

contains the costs for each point in relation to each possible 
track. The cost matrix is indexed along one axis by the track 
number, along another axis by the image number and along the 
third axis by a point number. Thus, each position in the cost 

10 matrix lists the cost for a unique combination of points and 
a track, one point from each image. Figure 2 also illustrates 
a {O, l} assignment matrix, z, which is defined with the same 
dimensions as the cost matrix. Setting a position in the 
assignment matrix to "one" means that the equivalent position 

15 in the cost matrix is selected into the solution. The 
illustrated solution matrix selects the {^24, 1122] solution 
previously described. Note that for the above example of two 
tracks and two snapshots, the resulting cost and assignment 
matrices are three-dimensional . As used in this patent 

20 application, the term "dimension" means the number of axes in 
the cost or assignment matrix while size refers to the number 
of elements along a typical axis. The costs and assignments 
have been grouped in matrices to facilitate computation. 
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A solution to the assignment problem satisfies two 
constraints - first, the sum of the associated costs for 
assigning points to a track extension is minimized and, 
second, if no false positives or false negatives exist, then 
5 each point is assigned to one and only one track. 

When false positives exist, however, additional 
hypothetical track extensions incorporating the false 
positives will be generated. Further note that the random 
locations of false positives will, in general, not fit well 

10 with true data and such additional hypothetical track 
extensions will result in higher costs. Also note that when 
false negative errors exist, then the size of the cost matrix 
must grow to include hypothetical track extensions formulated 
with "gaps" (i.e., data omissions where there should be 

15 legitimate observation data) for the false negatives. Thus, 
the second criteria must be weakened to reflect false 
positives not being as signed and also to permit the gap 
filler to be multiply assigned. With hypothetical cost 
calculated in this manner then the foregoing criteria for 

20 minimization will tend to materialize the false negatives and 
avoid the false positives. 

For a three-dimensional problem, as is illustrated in 
Figure 1, but with (initial) tracks, N2 observations in scan 
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1, N3 observations in scan 2, false positives and negatives 
assumed, the assignment problem can be formulated as: 

w, 

(a) Minimize: E E E ^1,^3^1,^3 

i^=0 12=0 13=0 123 

(b) Subject to: z2 z2 ^iii=^f i, = 1, . . .N,, (0.1) 
5 (c) z2 Z2 ^iii^^' i2 = I,...N2, 

i,=l i3=l ^ 2 3 

(d) £ 2-/ ^ili^l^ ^3 = 

i, = l i2 = l ^ 2 3 

{f.\ Z. . , E {0,1} V Z, , , , 



10 where "c" is the cost and "z" is a point or observation 
assignment, as in Figure 2. 

The minimization equation or equivalently objective 
function (0.1) (a) specifies the sum of the element by element 
product of the c and z matrices. The summation includes 

15 hypothesis representations z^^^^^^ with observation number zero 
being the gap filler observation. Equation (0.1) (b) requires 
that each of the tracks T^^.-.^T^^ be extended by one and only 
one hypothesis. Equation (0.1) (c) relates to each point or 
observation in the first observation set and requires that 

20 each such observation, except the gap filler, can only 
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associate with one track but, because of the "less than" 
condition, it might not associate with any track. Equation 
(0.1) (d) is like (0.1) (c) except that it is applicable to the 
second observation set. Equation (0.1) (e) requires that 
elements of the solution matrix z be limited to the zero and 
one values. 

The only known method to solve Problem Formulation (0.1) 
exactly is a method called "Branch and Bound." This method 
provides a systematic ordering of the potential solutions so 
that solutions with a same partial solution are accessible via 
a branch of a tree describing all possible solutions whereby 
the cost of unexamined solutions on a branch are incrementally 
developed as the cost for other solutions on the branch are 
determined. When the developing cost grows to exceed the 
previously known minimal cost (i.e., the bound) then 
enumeration of the tree branch terminates. Evaluation 
continues with a new branch. If evaluation of the cost of a 
particular branch completes, then that branch has lower cost 
than the previous bound so the new cost replaces the old 
bound. When all possible branches are evaluated or eliminated 
then the branch that had resulted in the last used bound is 
the solution. If we assume that = N2 = ^ n and that 
branches typically evaluate to half there full length, then 
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workload associated with "branch and bound" is proportional to 
(nil — !)^. This workload is unsuited to real time evaluation. 

The Branch and Bound algorithm has been used in past 
research on the Traveling Salesman Problem. Messrs. Held and 
5 Karp showed that if the set of constraints was relaxed by a 
method of Lagrangian Multipliers (described in more detail 
below) then tight lower bounds could be developed in advance 
of enumerating any branch of the potential solution. By 
initiating the branch and bound algorithm with such a tight 
10 lower bound, significant performance improvements result in 
that branches will typically evaluate to less than half their 
full length. 

Messrs. Frieze and Yadagar in dealing with a problem 
related to scheduling combinations of resources, as in job, 

15 worker and work site, showed that Problem Formulation (0.1) 
applied. They further described a solution method based upon 
an extension of the Lagrangian Relaxation previously 
mentioned. The two critical extensions provided by Messrs. 
Frieze and Yadagar were: (1) an iterative procedure that 

2 0 permitted the lower bound on the solution to be improved (by 
"hill climbing" described below) and (2) the recognition that 
when the lower bound of the relaxed problem was maximized, 
then there existed a method to recover the solution of the 
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non-relaxed problem in most cases using parameters determined 
at the maximum. The procedures attributed to Messrs. Frieze 
and Yadagar are only applicable to the three-dimensional 
problem posed by Problem Formulation (0.1) and where the cost 



airborne objects usually requires solution of a much higher 
dimensional problem. 

Figures 1 and 2 illustrate an example where "look ahead" 
data from the second image improved the assignment accuracy 

10 for the first image. Without the look ahead, and based only 
upon a simple nearest neighbor approach, the assignments in 
the first set would have been reversed. Problem Formulation 
(0.1) and the prior art only permit looking ahead one image. 
In the prior art it was known that the accuracy of assignments 

15 will improve if the process looks further ahead, however no 
practical method to optimally incorporate look ahead data 
existed. Many real radar tracking problems involve hundreds 
of tracks, thousands of observations per observation set and 
matrices with dimensions in the range of 3 to 2 5 including 

2 0 many images of look ahead. 

It was also known that the data assignment problem may be 
simplified (without reducing the dimension of the assignment 
problem) by eliminating from consideration for each track 



5 matrix is fully populated. 



However, tracking multiple 
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those points which, after considering estimated limits of 
speed and turning ability of the objects, could not physically 
be part of the track. One such technique, denoted hereinafter 
the "cone method, " defines a cone as a continuation of each 
5 previously determined track with the apex of the cone at the 
end of the previously defined track. The length of the cone 
is based on the estimated maximum speed of the object and the 
size of the arc of the cone is based on the estimated maximum 
turning ability of the object. Thus, the cone defines a 

10 region outside of which no point could physically be part of 
the respective track. For any such points outside of the 
cones, an infinite number could be put in the cost matrix and 
a zero could be preassigned in the assignment matrix. It was 
known for the tracking problem that these elements will be 

15 very common in the cost and selection matrices (so these 
matrices are "sparse"). 

It was also known in the prior art that one or more 
tracks which are substantially separated geographically from 
the other tracks can be separated also in the assignment 

20 problem. This is done by examining the distances from each 
point to the various possible tracks. If the distances from 
one set of points are reasonably short only in relation to one 
track, then they are assigned to that track and not further 
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considered with the remainder of the points. Similarly, if a 
larger group of points can only be assigned to a few tracks, 
then the group is considered a different assignment problem. 
Because the complexity of assignment problems increases 
5 dramatically with the number of possible tracks and the total 
number of points in each matrix, this partitioning of the 
group of points into a separate assignment problem and removal 
of these points from the matrices for the remaining points, 
substantially reduces the complexity of the overall assignment 
10 problem. 

A previously known Multiple Hypothesis Testing (MHT) 
algorithm (see Chapter 10 of S.S, Blackman. Multiple-Target 
Tracking with Radar Applications. Artech House, Norwood, MA, 
1986.) related to formulation and scoring. The MHT procedure 

15 describes how to formulate the sparse set of all reasonable 
extension hypothesis (for Figure 1 the set {h^2---^24}) how 
to calculate a cost of the hypothesis [T^, o^j, 02j^} based upon 
the previously calculated cost for hypothesis {T^, Oj^} . The 
experience with the MHT algorithm, known in the prior art, is 

2 0 the basis for the assertion that look ahead through k sets of 
observations results in improved assignment of observations 
from the first set to the track. 
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In theory, the MHT procedure uses the extendable costing 
procedure to defer assignment decision until the accumulated 
•evidence supporting the assignment becomes overwhelming. When 
it makes the assignment decision, it then eliminates all 
5 potential assignments invalidated by the decision in a process 
called "pruning the tree." In practice, the MHT hypothesis 
tree is limited to a fixed number of generations and the 
overwhelming evidence rule is replaced by a most likely and 
feasible rule. This rule considers each track independently 
10 of others and is therefore a local decision rule. 

A general object of the present invention is to provide 
an efficient and accurate process for assigning each point 
object in a region from multiple images to a proper track and 
then taking an action based upon the assignments. 
15 A more specific object of the present invention is to 

provide a technique of the foregoing type which determines the 
solution of a ^c-dimensional assignment problem where "k" is 
greater than or equal to three. 

2 0 SUMMARY OF THE INVENTION 

The present invention relates to a method and apparatus 
for tracking objects. In particular, the present invention 
tracks movement or trajectories of objects by analyzing 
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radiation reflected from the objects, the invention being 
especially useful for real-time tracking in noisy 
environments . 

In providing such a tracking capability, a region 
5 containing the objects is repeatedly scanned to generate a 
multiplicity of sequential images or data observation sets of 
the region. One or more points (or equivalent ly 

observations) , in each of the images or observation sets are 
detected wherein each such observation either corresponds to 

10 an actual location of an object or is an erroneous data point 
due to noise. Subsequently, for each observation detected, 
figures of merit or costs are determined for assigning the 
observation to each of a plurality of previously determined 
tracks. *Afterwards, a first optimization problem is 

15 specified which includes: 

(a) a first objective function for relating the above 
mentioned costs to potential track extensions through the 
detected observations (or simply observations) ; and 

(b) a first collection of constraint sets wherein each 
20 constraint set includes constraints to be satisfied by 

the observations in a particular scan to which the 
constraint set is related. In general, there is a 
constraint for each observation of the scan, wherein the 
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constraint indicates the number of track extensions to 
which the observation may belong. 

In particular, the first optimization problem is 
formulated, generated or defined as an M-dimensional 
5 assignment problem wherein there are M constraint sets in the 
first collection of constraint sets (i.e., there are M scans 
being examined) and the first objective function minimizes a 
total cost for assigning observations to various track 
extensions wherein terms are included in the cost, such that 

10 the terms have the figures of merit or costs for hypothesized 
combinations of assignments of the observations to the 
tracks. Subsequently, the formulated M-dimensional assignment 
problem is solved by reducing the complexity of the problem by 
generating one or more optimization problems each having a 

15 lower dimension and then solving each lower dimension 
optimization problem. That is, the M-dimensional assignment 
problem is solved by solving a plurality of optimization 
problems each having a lower number of constraint sets. 

The reduction of the M-dimensional assignment problem to 

20 a lower dimensioned problem is accomplished by relaxing the 
constraints on the points of one or more scans thereby 
permitting these points to be assigned to more than one track 
extension. In relaxing the constraints, terms having penalty 
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factors are added into the objective function thereby 
increasing the total cost of an assignment when one or more 
points are assigned to more than one track. Thus, the 
reduction in complexity by this relaxation process is 
5 iteratively repeated until a sufficiently low dimension is 
attained su'ch that the lower dimensional problem may be solved 
directly by known techniques. 

In one embodiment of the invention, each Jc-dimensional 
assignment problem [2 < k < M) is iteratively reduced to a {k~ 

10 1) - dimensional problem until a two-dimensional problem is 
specified or formulated. Subsequently, the two-dimensional 
problem formulated is solved directly and a "recovery" 
technique is used to iteratively recover an optimal or near- 
optimal solution to each A:-dimensional problem from a derived 

15 (A:-l) -dimensional problem k= 2,3,4,...M. 

In performing each recovery step (of obtaining a solution 
to a Jc-dimensional problem using a solution to a {k-1) - 
dimensional problem) , an auxiliary function is utilized. In 
particular, to recover an optimal or near-optimal solution to 

20 a k-dimensional problem, an auxiliary function, W^^^^ 
CA:=4,5, . . . ,M) , is specified and a region or domain is 
determined wherein this function is maximized, whereby values 
of the region determine the penalty factors of the {k-1) - 
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dimensional problem such that another two-dimensional problem 
can be formulated which determines a solution to the k- 
dimensional problem using the penalty factors of the {k~l) - 
dimensional problem . 
5 Each Auxiliary function iPj,, ;c=3 , . . . , M-1 , is a function 

of both lower dimensional problems, penalty factors, and a 
solution at the dimension k at which the penalized cost 
function is solved directly (typically a two-dimensional 
problem) . Further, in determining, for auxiliary function W^, 

10 a peak region, a gradient of the auxiliary function is 
determined, and utilized to identify the peak region. Thus, 
gradients are used for each of the approximation functions W^, 
W^, , . . , ¥^,1 which are sequentially formulated and used in 
determining penalty factors until M-1 is used in determining 

15 the penalty factors for the (M-1) -dimensional problem. 
Subsequently, once the M-dimensional problem is solved (using 
a two-dimensional problem to go from an (M-1) -dimensional 
solution to an M-dimensional solution) , one or more of the 
following actions are taken based on the track assignments: 

20 sending a warning to aircraft or a ground or sea facility, 
controlling air traffic, controlling anti-aircraft or anti- 
missile equipment, taking evasive action, working on one of 
the objects. 
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According to one feature of this first embodiment of the 
present invention, the following steps are also performed 
before the step of defining the auxiliary function. A 
preliminary auxiliary function is defined for each of the 
5 lower dimensional problems having a dimension equal or one 
greater than the dimension at which the penalized cost 
function is solved directly. The preliminary auxiliary 
function is a function of lower order penalty values and a 
solution at the dimension at which the penalized cost function 

10 was solved directly. In determining a gradient of the 
preliminary auxiliary function, step in the direction of the 
gradient to identify a peak region of the preliminary 
auxiliary function and determine penalty factors at the peak 
region. Iteratively repeat the defining, gradient 

15 determining, stepping and peak determining steps to define 
auxiliary functions at successively higher dimensions until 
the auxiliary function dimension {k-1) is determined. In an 
alternative second embodiment of the present invention, 
instead of reducing the dimentionality of the M-dimensional 

20 assignment problem by a single dimension at a time, a 
plurality of dimensions are relaxed simultaneously. This new 
strategy has the advantage that when the M-dimensional problem 
is relaxed directly to a two-dimensional assignment problem, 



18 



CSUR.01USRI 

then all computations may be performed precisely without 
utilizing an auxiliary function such as as in the first 
embodiment. More particularly, the second embodiment solves 
the first optimization problem (i.e., the M-dimensional 
assignment problem) by specifying (i.e., creating, generating, 
formulating and/or defining) a second optimization problem. 
The second optimization problem includes a second objective 
function and a second collection of constraint sets wherein: 

a) the second objective function is a combination of the 
first objective function and penalty factors or terms 
determined for incorporating {M-m) -constraint sets of the 
first optimization problem into the second objective 
function; 

b) the constraint sets of the second collection include only 
m of the constraint sets of the first collection of 
constraints, 

wherein 2 m < M ~1 . Note that, once the second optimization 
problem has been specified or formulated, an optimal or near- 
optimal solution is determined and that solution is used in 
specifying (i.e, creating, generating, formulating and/or 
defining) a third optimization problem of {M-m) dimensions (or 
equivalently constraint sets) . The third optimization problem 
is subsequently solved by decomposing it using the same 
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procedure of this second embodiment as was used to decompose 
the first optimization problem above. Thus, a plurality of 
instantiations of the third optimization problem are 
specified, each successive instantiation having a lower number 
5 of dimensions, until an instance of the third optimization 
problem is a two-dimensional assignment problem which can be 
solved directly. Subsequently, whenever an instance of the 
third optimization problem is solved, the solution is used to 
recover a solution to the instance of the first optimization 

10 problem from which this instance of the third optimization was 
derived. Thus, an optimal or near-optimal solution to the 
original first optimization problem may be recovered through 
iteration of the above steps. 

As mentioned above, the second embodiment of the present 

15 invention is especially advantageous when m =2, since in this 
case all computations may be performed precisely and without 
utilizing auxiliary functions. 

Other features and benefits of the present invention will 
become apparent from the detailed description with the 

2 0 accompanying drawings contained hereinafter. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 is a graph of images or data sets generated by 
a scan of a region and possible tracks within the images or 
data sets according to the prior art. 

Figure 2 illustrates cost and assignment matrices for the 
data sets of Figure 1 according to the prior art. 

Figure 3 is a block diagram of the present invention. 

Figure 4 is a flow chart of a process according to the 
prior art for solving a three-dimensional assignment problem. 

Figure 5 is a flow chart of a process according to the 
present invention for solving a A:-dimensional assignment 
problem where "A:" is greater than or equal to 3. 

Figure 6 is a graph of various functions used to explain 
the present invention. 

Figure 7 is another block diagram of the present 
invention for solving the Jc-dimensional assignment problem 
where "k" is greater than or equal to three (3) . 

Figure 8 is a flowchart describing the procedure for 
solving an n-dimensional assignment problem according to the 
second embodiment of the invention. 
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DETAILED DESCRIPTION OF THE INVENTION 
Referring now to the other figures in detail wherein like 
reference numerals indicate like elements throughout the 
several views, Figure 3 illustrates a system generally 
5 designated 100 for implementing the present invention. System 
100 comprises, for example, a radar station 102 (note sonar, 
microwave, infrared and other radiation bandwidths are also 
contemplated) for scanning a region which may be, for example, 
an aerial region (in aerial surveillance applications) or a 

10 work region (in robotic applications) and generating signals 
indicating locations of objects within the region. The 
signals are input to a converter 104 which converts the 
signals to data points or observations in which each object 
(or false positive) is represented by a single point. The 

15 output of the converter is input to and readable by a computer 
106. As described in more detail below, the computer 106 
assigns the points to respective tracks, and then displays the 
tracks and extrapolations of the tracks on a monitor 110. 
Also, the computer 106 determines an appropriate action to 

20 take based on the tracks and track extensions. For example, 
in a commercial application at an airport, the computer can 
determine if two aircraft being tracked are on a collision 
course and if so, signal a transmitter 112 to warn each 
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aircraft, or if a scheduled take-off will pose the risk of 
collision, delay the take-off. For a military application on 
a ship or base, the computer can determine subsequent 
coordinates of enemy aircraft and send the coordinates to an 
5 antiaircraft gun or missile 120 via a communication channel 
122. In a robotic application, the computer controls the 
robot to work on the proper object or portion of the object. 
The invention generates Jc-dimensional matrices where k is 
:;f the number of images or sets of observation data in the look 

y 

J;;; 10 ahead window plus one. Then, the invention formulates a k- 
p dimensional assignment problem as in (0.1) above, 

p The A:-dimensional assignment problem is subsequently 

relaxed to a {k-1) -dimensional problem by incorporating one 

•AXi 

\ \ 

[j set of constraints into the objective function using a 

15 Lagrangian relaxation of this set. Given a solution of the 
(^-1) -dimensional problem, a feasible solution of the k- 
dimensional problem is then reconstructed. The {k-1) - 
dimensional problem is solved in a similar manner, and the 
process is repeated until it reaches the two-dimensional 
2 0 problem that can be solved exactly. The ideas behind the 
Lagrangian relaxation scheme are outlined next. 
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Consider the integer programming problem 
Minimize v(z) = c^z 

Subject to: Az < h, (0.2) 
5 Bz < d, 

Zi is an integer for i E I, 

where the partitioning of the constraints is natural in some 

sense. Given a multiplier vector u > 0, the Lagrangian 

10 relaxation of (0.2) relative to the constraints Bz < d is 

defined to be 

0(u) = Minimize (c'^z-^u^ (Bz-d) } (0.3) 

Subject to: Az < b, 
15 is an integer for i € I. 

If the constraint set Bz < d is replaced by Bz = d, the 

non-negativity constraint on u is removed. £f=c'^z-hu'^ (Bz-d) is 

the Lagrangian relative to the constraints Bz < d, and hence 
2 0 the name Lagrangian relaxation. The following fact gives the 

relationship between the objective functions of the original 

and relaxed problems. 

FACT A.l. If ¥ is an optimal solution to (0.2), then 

0(u)^ vCz) for all u ^ 0. If an optimal solution z of (0.3) 
25 is feasible for (0.3),. then z is an optimal solution for (0.2) 

and <P(u) = v( z) . 

This leads to the following algorithm: 
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Algorithm. Construct a sequence of multipliers {u^]k=o 
converging to the solution u" of Maximize {^(u) : u > O} and a 
corresponding sequence of feasible solutions {z^]'^^q of (0.2) as 
follows : 

5 1. Generate an initial approximation u^. 

2. Given u^^, choose a search direction s^^ and a search 
distance so that <P{Ui, ^a^^s^) > 0{uj,) . Update the multiplier 
estimate U;^^ by Uj,^-^ = U;, + ctj^Sj^, 

n 3. Given Uj,^^ and a feasible solution z^^iiuj^^i) of (0.3), 

y 

fi 10 recover a feasible solution Zj,^i{uj,^i) of the integer programming 

? 

% problem (0.2) . 

h 

y 4. Check the termination criteria. If the algorithm is 

'I not finished, set k=k^l and return to Step 2. Otherwise, 

U terminate the algorithm. 

:i3 15 If z" is an optimal solution of (0.2), then 

0(U},) < 0(u) ^ viz) < v(zj . 
Since the optimal solution z = z(u) of (0.3) is usually not a 
feasible solution of (0.2) for any choice of the multipliers, 
0(u) is usually strictly less than vfz). The difference vfz) - 
20 0(u) is called the "duality gap, " which can be estimated by 
comparing the best relaxed and recovered solutions computed in 
the course of maximizing 0(u) . 
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Thus, to eliminate the "<" in the constraint formulations 
(e.g., (0.1) (a) - (0.1) (d)) that resulted from false 
positives, the constraint sets have been modified, as 
indicated above, to incorporate a gap filler in each 
5 observation set to account for false positives. Thus, a zero 
for an index {1 < p < k) ±n in problem (0.4) below 

indicates that the track extension represented by ^i^^...!^ 
includes a false positive in the p^^ observation set. Note 
that this implies that a hypothesis be formed incorporating an 
10 observation with k-1 gap fillers, e.g., ZQ^^Q^^o^^ot ip^O- Thus, 



the resulting generalization of Problem Formulation (0.1) 



the following Jc-Dimensional Assignment Problems in which k >3: 



without the "less than" complication within the constraints is 
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Minimize E • • • i^i i 



for ij = 1, . . .Nj 
and j=2 , . . .k - 1, 



z. . e {O, l}, W-i n-7 Jk- 



10 where c and z are similarly dimensioned matrices representing 
costs and hypothetical assignments. Note, in general, for 
tracking problems these matrices are sparse. 

After formulating the Jc-dimensional assignment problem as 
in (0.4), the present invention solves the resulting problem 

15 so as to generate the outputs required by devices 110, 112, 
122 and 130. For each observation set 0^ received from 
converter 104 at time where 1=1,.,,, N^, the computer 

processes in a batch together with the other observation 
sets Oi_j^^i, . . . , Oi and the track Ti.^, i.e., Tj is the set of all 

20 tracks that have been defined up to but not including Oj . 
(Note, bold type designations refer to the vector of elements 
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for the indicated time, i.e., the set of all observations in 
the scan or tracks existing at the time, etc.) The result of 
this processing is the new set of tracks T^.j^^i and a set of 
cost weighted possible solutions indicating how the tracks 
5 might extend to the current time , At time t^^^ the batch 
process is repeated using the newest observation set and 
deleting the oldest. Thus, there is a moving window of 
observation sets which is shifted forward to always 
incorporate the most recent observation set. The effect is 

10 that input observation sets are reused for (k-1) cycles and 
then on the observation set ^ s k^^ reuse each observation within 
the observation set is integrated into a track. 

Figure 7 illustrates various processes implemented upon 
receipt of each observation set. Except for the addition of 

15 the Tc-dimensional assignment solving process 3 00 and the 
modification to scoring process 154 to build data structures 
suitable for process 300, all processes in Figure 7 are based 
upon prior art. The following processes 150 and 152 extend 
previously defined tracks h^.i based on new observations. Gate 

20 formulation and output process 156 determines, for each of the 
previously defined tracks, a zone wherein the track may 
potentially extend based on limits of velocity, 
maneuverability and radar precision. One such technique to 
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accomplish this is the cone method described previously. The 
definition of the zone is passed to gating process 150. When 
a new observation set Oi is received, the gating process 150 
will match each member observation with the zone for each 
5 member of the hypothetical set h^.i. After all input 
observations from Oi are processed, the new hypothesis set 
is generated by extending each track of the prior set of 
hypothetical tracks hi.^ either with missed detect gap fillers 
or with all new observation elements satisfying the track's 

10 zone. This is a many-to-many matching in that each hypothesis 
member can be extended to many new observations and each new 
observation can be used to extend many hypotheses. It, 
however, is not a full matching in that any hypothesis will 
neither be matched to all observations nor vice versa. It is 

15 this matching characteristic that leads to the sparse matrices 
involved in the tracking process. Subsequently, gating 150 
forwards the new hypothesis set hi to filtering process 152. 
Filtering process 152 determines a smooth curve for each 
member of hi. Such a smooth curve is more likely than a sharp 

2 0 turn from each point straight to the next point. Further, the 
filtering process 152 removes small errors that may occur in 
generating observations. Note that in performing these tasks, 
the filtering process 152 preferably utilizes a minimization 

29 
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of a least squares test of the points in a track hypothesis or 
a Kalman filtering approach. 

As noted above, the foregoing track extension process 
requires knowledge of a previous track. For the initial 
5 observations, the following gating process 158 and filtering 
process 160 determine the "previous track" based on initial 
observations. In determining the initial tracks, the points 
from the first observation set form the beginning points of 
all possible tracks. After observation data from the next 

10 observation set is received, sets of simple two-point straight 
line tracks are defined. Then, promotion, gate formulation, 
and output step 162 determines a zone in which future 
extensions are possible. Note that filtering step 160 uses 
curve fitting techniques to smooth the track extensions 

15 depending upon the number of prior observations that have 
accumulated in each hypothesis. Further note that promotion, 
gate formulation and output process 162 also determines when 
sufficient observations have accumulated to form a basis for 
promoting the track to processes 150 and 152 as described 

2 0 above . 

The output of each of the filtering processes 152 and 160 
is a set of hypothetical track extensions. Each such 
extension contains the hypothetical initial conditions (from 
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the previous track) , the list of observations incorporated in 
the extension, and distance between each observation and the 
filtered track curve. Scoring process 154 determines the 
figure of merit or cost of an observation being an extension 
5 of a track. In one embodiment, the cost is based on the 
above-mentioned distance although the particular formula for 
determining the cost is not critical to the present invention. 
A preferred formula for determining the cost utilizes a 
negative log likelihood function in which the cost is the 

10 negative of the sum of: (a) the logs of the distances 
normalized by sensor standard deviation parameters, and (b) 
the log likelihoods for events related to track initiation, 
track termination, track maneuver, false negatives and false 
positives- Note that track maneuvers are detected by 

15 comparing the previous track curve with the current extension. 
Further note that some of the other events related to, for 
example, false negatives and false positives are detected by 
analyzing the relative relationship of gap fillers in the 
hypothesis. Thus, after determining that one of these events 

20 occurred, a cost for it can be determined based upon suitable 
statistics tables and system input parameters. The negative 
log likelihood function is desirable because it permits 
effective integration of the useful components. Copies of the 
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set of hypothetical track extensions which are scored are 
subsequently passed directly to one of the gate formulation 
and output steps 156 and 162. Note that the scoring process 
154 also arranges the actual scores in a sparse matrix based 
5 upon observation identifiers, and passes them to A:-dimensional 
assignment problem solving process 300. 

The assignment solving process 300 is described below. 
Its output is simply the list of assignments which constitute 
the most likely solution of the problem described by Equation 

10 (0.4). Note that both gate formulation and output processes 
156 and 162 use (at different times) the list of assignments 
to generate the updated track history Ti to eliminate or prune 
alternative previous hypotheses that are prohibited by the 
actual assignments in the list, and, subsequently, to output 

15 any required data. Also note that when one of the gate 
formulation and output processes 156 and 162 accomplish these 
tasks, the process will subsequently generate and forward the 
new set of gates for each remaining hypothesis and the 
processes will then be prepared to receive the next set of 

2 0 observations. In one embodiment, the loop described here will 
generate zones for a delayed set of observations rather than 
the subsequent set. This permits processes 156 and 162 to 
operate on even observation sets while the scoring step 154 
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and A:-dimensional solving process operate on odd sets of 
observations, or vice versa. 

The assignment solving process 300 permits the present 
invention to operate with a window size of dimension k-1 for 
5 some k > 3. The upper limit on k depends only upon the 
computational power of the computer 106 and the response time 
constraints of system 100. The ^c-1 observation sets within 
the processing window plus the prior track history result in 
a ^c-dimensional Assignment Problem as described by Problem 
10 Formulation (0.4). The present invention solves this 
generalized problem including the processes required to 
consider false positives and negatives, and also the processes 
required to consider sparse matrix problem formulations. 

15 I, A FIRST EMBODIMENT OF THE 

;c-DIMENSIONAL ASSIGNMENT SOLVER 300 

In describing a first embodiment of the Jc-dimensional 

assignment solver 300, it is worthwhile to also discuss the 

process of Figure 4 which is used by the solver 300. Figure 

2 0 4 illustrates use of the Frieze and Yadagar process as shown 

in prior art for transforming a three-dimensional assignment 

problem into a two-dimensional assignment problem and then use 

a hill climbing algorithm to solve the three-dimensional 

assignment problem. The solver 3 00 uses a Lagrangian 
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Relaxation technique (well known in the art) to reduce the 
dimension of an original ^-dimensional assignment problem 

{k > 3) down to a three-dimensional problem and then use the 
process of Figure 4 to solve the three-dimensional problem. 
Further note that the Lagrangian Relaxation technique is also 
utilized by the process of Figure 4 and that in using this 
technique the requirement that each point is assigned to one 
and only one track is relaxed. Instead, an additional cost, 
which is equal to a respective Lagrangian Coefficient u, is 
added to the cost or objective function (0.1) (a) whenever a 
point is assigned to more than one track. This additional 
cost can be picked to weight the significance of each 
constraint violation differently, so this additional cost is 
represented as a vector of coefficients u which are correlated 
with respective observation points. Hill climbing will then 
develop a sequence of Lagrangian Coefficients sets designated 

(Uq, . . . ,Uj, Uj+i, . . . ,Up) . That correspond to an optimum solution 
of the two-dimensional assignment problem. The assignments at 
this optimum solution are then used to "recover" the 
assignment solution of the three-dimensional assignment 
problem. 

In step 200 of Figure 4, initial values are selected for 
the Uq coefficients. Because the Lagrangian Relaxation process 
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is iterative, the initial values are not critical and are all 
initially selected as zero. In step 202, these additional 
costs are applied to the objective function (0.1) (a). With 
the addition of the costs u, the goal is still to assign the 
points which minimize the total cost. This transforms 
Equation (0.1) (a), written for 3 and altered to exclude 
mechanisms related to false positives and negatives, into 
objective function (1.1) (a) . In the first iteration it is not 
necessary to consider the u matrix because all u values are 
set to zero. To relax the requirement that each point be 
assigned to one and only one track, the constraint Equation 
(0.1) (d) is deleted, thereby permitting points from the last 
image to be assigned to more than one track. Note that while 
any axis can be chosen for relaxation, observation constraints 
are preferably relaxed. The effect of this relaxation is to 
create a new problem which must have the same solution in the 
first two axes but which can have a differ solution in the 
third axis. The result is constraints (1.1) (b-d) . 
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(a) 



Minimize 




(b) Subject to: 




1, . . 



• / 



(1.1) 



,(c) 



(d) 




^1/ ±2 



V z 



'2f 



5 



Step 2 04 then generates from the three-dimensional problem 
described by Problem Formulation (0.1) a new two-dimensional 
problem formulation which will have the same solution for the 

10 first two indices. As Problem Formulation (1.1) has no 
constraints on the 3"^^ axis, any value within a particular 3"^^ 
axis can be used in a solution, but using anything other than 
the minimum value from any 3"^^ axis has the effect of 
increasing solution cost. Conceptually, the effect of step 

15 2 04 is to change the three-dimensional arrays ^ in Problem 
Formulation (1.1) into two-dimensional arrays as shown in 
Problem Formulation (1.2) and * to generate the new two- 
dimensional matrix ^i^i^ defined as shown in Equation (1.3) . 
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(a) Minimize 

(b) Subject to: 

(c) 
(d) 

5 

m. . = Minimum {c..-u. |2 = l,...,Ar. }. (1.3) 

The cost or objective function for the reduced problem as 
defined by (1.2) (a), if evaluated at all possible values of u 
is a surface over the domain of u. This surface is referred 

10 to as ^(u) and is non-smooth but provable convex (i.e., it has 
a single peak and several other critical characteristics which 
form terraces) . Due to the convex characteristics of ^(u) , 
the results from solving Problem Formulation (1.2) at any 
particular Uj can be used to generate a new set of coefficients 

15 Uj^i whose corresponding cost value is closer to the peak of 
^(u) . The particular set of Lagrangian Coefficients that will 
generate the two-dimensional problem resulting in the maximum 
cost is designated vip. To recover the solution to the three- 
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i,=i 23=1 2 3 



= l i3 = l 

z., e {0,1} V z, 



CSUR.01USRI 

dimensional assignment problem requires solving the Equation 

(1.2) problem corresponding to Up. 

In step 2 06, the two-dimensional problem is solved 
directly using a technique known to those skilled in the art 
such as Reverse Auction for the corresponding cost and 
solution values. This is the evaluation of one point on the 
surface or for the first iteration ^'(Uo) - 

Thus, after this first iteration, the points have been 
assigned based on all "u" values being arbitrarily set to 
zero. Because the "u" values have been arbitrarily assigned, 
it is unlikely that these assignments are correct and it is 
likely that further iterations are required to properly assign, 
the points. Step 208 determines whether the points have been 
properly assigned after the first iteration by determining if 
for this set of assignments whether a different set of "u" 
values could result in a higher total cost. Thus, step 208 is 
implemented by determining the gradient of objective function 

(1.2) (a) with respect to Uj . If the gradient is substantially 
non-zero (greater than a predetermined limit) then the 
assignments are not at or near enough to the peak of the 
surface ^'(u) (decision 210) , and the new set of Lagrangian 
Coefficients Uj^^ is determined. 
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Hill climbing Step 212 determines the Uj^^ values based 
upon the Uj values, the direction resulting from protecting the 
previous gradient into the domain, and a step size. The 
solution value of the two-dimensional problem is the set of 
5 coefficients that minimize the two-dimensional problem and the 
actual cost at the minimum. Those coefficients augmented by 
the coefficients stored in m^^^^ permit the evaluation (but not 
the minimization) of the cost term in Problem Formulation 

(1.1) . These two cost terms are lower and upper bounds on the 
10 actual minimized cost of the three-dimensional problem, and 

the difference between them in combination with the gradient 
is used to compute the step size. 

With this new set of u values, steps 202-210 are repeated 
as a second iteration. Steps 212 and 202-210 are repeated 
15 until the gradient as a function of u determined in step 2 08 
is less than the predetermined limit. This indicates that the 
Tip values which locate the peak area of the ^(u) surface are 
determined and that the corresponding Problem Formulation 

(1.2) has been solved. Step 214 will attempt to use the 
20 assignments that resulted from this particular two-dimensional 

assignment problem to recover the solution of the three- 
dimensional assignment problem as described below. If the 
limit was chosen properly so that the u values are close 
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enough to the peak, this recovery will yield the proper set of 
assignments that rigidly satisfies the constraint that each 
point be assigned to one and only one track. However, if the 
u values are not close enough to the peak, then the limit 
value for decision 210 is reduced and the repetition of steps 
212 and 202-210 is continued. 

Step 214 recovers the three-dimensional assignment 
solution by using the assignment values determined on the last 
iteration through step 208. Consider the two-dimensional z 
assignment matrix to have I's in the locations specified by 
the list L2=(ai, bi)^^^. If the three-dimensional z matrix is 
specified by placing 1 * s at the location indicated by the list 
Lp=(ai, h^, inaiJbi)i=i then the result is a solution of Problem 
Formulation (1.1). Let 1/3= (/n^^jb^) ^^-^ be the list formed by the 
third index. If each member of L3 is unique then the L2 
solution satisfies the third constraint so it is a solution to 
Problem Formulation (0.1). When this is not the case, 
recovery determines the minimal substitutions required within 
list L3 so that it plus L^^ will be a feasible solution, i.e., 
a solution which satisfies the constraints of a problem 
formulation, but which may not optimize the objective function 
of the problem formulation. This stage of the recovery 
process is formulated as a two-dimensional assignment problem 
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as follows. Form a new cost matrix [Ci^^] i^^^^ where Ci^^^c^^^,^^^ for 
j=l...i\^i and the term is the total number of cost elements 
in the selected row of the three-dimensional cost matrix. 
Attempt to solve this two-dimensional problem for the minimum 
5 using two constraints sets. If a feasible solution is found 
then the result will have the same form as list L^. Replace 
the first set of indices by the indicated (a^, h^) pairs taken 
from list and the result will be a feasible solution of 
Problem Formulation (0.4) . If no feasible solution to the new 
10 two-dimensional problem exists then further effort to locate 
the peak of ^(u) is required . 



I . 1 . Generalization to a Multi-Dimensional 
Assignment SolvinQ Process 

15 Let M be a fixed integer and assume that M is the 

dimension of the initial assignment problem to be solved. 
Thus, initially, the result of the scoring step 154 is a M- 
dimensional Cost Matrix which is structured as a sparse matrix 
(i.e., only a small percentage of the entries in the cost and 

20 assignment matrices are filled or non-zero) . Individual cost 
elements represent the likelihood that a track as extended 
by the set of observations {0^^-! i = l,...,M-l}, is not valid. 
Because the matrix is sparse the list of cost elements is 
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stored as a packed list, and then for each dimension of the 
matrix, a vector of a variable length list of pointers to the 
cost elements is generated and stored. This organization 
means that for a particular observation O^j for the j^^ list in 
5 the i^^ vector will be a list of pointers to all hypotheses in 
which Oij- participates. This structure is further explained in 
the following section dealing with problem partitioning. 

The objective of the assignment solving process is to 
select from the set of all possible combinations of track 

10 extensions a subset that satisfies two criteria. First, each 
point in the subset of combinations should be assigned to one 
and only one track and therefore, included in one and only one 
combination of the subset, and second, the total of the 
scoring sums for the combinations of the subset should be 

15 minimized. This yields the following M-dimensional equations 
where, k = M: 
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(a) Minimize v, 




(1.4) 



5 



(b) Subject to: 




= I 



/ • • • 



'1/ 



(c) 




i,=0 i,.,=0 i,,,=0 i,=0 



(d) 




f * • • t ^ y^i 



(e) 



z]" ,6(0,1} for all i 



1' * ' • ' -^k^ 



and where c'^ is the cost matrix [c^^...i^] which is a function of 
15 the distance between the observed point and the smoothed 
track determined by the filtering step, and is the cost 



presented in Problem Formulation (0.4) except that it includes 
the subscript and superscript k notation. Thus, in solving 

2 0 the AT- dimensional Assignment Problem the invention reduces 
this problem to an (M-1) -dimensional Assignment Problem and 
then to an (i\r-2) -dimensional Assignment Problem, etc. 
Further, the symbol A:e{3, . . . ,m} customizes Problem Formulation 
(1.4) to a particular relaxation level. That is, the notation 

2 5 is used to reference data from levels relatively removed as in 
^k+i cost coefficients which existed prior to this level 



function. 



This set of equations is similar to the set 
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of relaxed coefficients c^. Note that actual observations are 
numbered from 1 to N^, where is the number of observations 
in observation set i. Further note that the added O 
observation in each set of observations is the unconstrained 
"gap filler." This element serves as a filler in substituting 
for missed detects, and a sequence of these elements including 
only one true observation represents the possibility that the 
observation is a false positive. Also note that by being 
unconstrained a gap filler may be used in as many hypotheses 
as required. 

While direct solution to (1.4) would give the precise 
assignment, the solution of Jc-dimensional equations directly 
for large k is too complex and time consuming for practice. 
Thus, the present invention solves this problem indirectly. 

The following is a short description of many aspects of 
the present invention and includes some steps according to the 
prior art. The first step in solving the problem indirectly 
is to reduce the complexity of the problem by the previously 
known and discussed Lagrangian Relaxation technique. 
According to the Lagrangian Relaxation technique, the absolute 
requirement that each point is assigned to one and only one 
track is relaxed such that for some one image, points can be 
assigned to more than one track. However, a penalty based on 
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a respective Lagrangian Coefficient is added to the cost 
function when a point in the image is assigned to more than 
one track. The Lagrangian Relaxation technique reduces the 
complexity or "dimension" of the formulation of the assignment 
5 problem because constraints on one observation set are 
relaxed. Thus, the Lagrangian Relaxation is performed 
iteratively to repeatedly reduce the dimension until a two- 
dimensional penalized cost function problem results as in 
Problem Formulation (1.1). This two-dimensional problem is 

n 

10 solved then directly by a previously known technique such as 

hi 
'^'\ 

^ Reverse Auction, The penalized cost function for the two- 

^« dimensional problem defines a valley or convex shaped surface 

which is a function of various sets of {u'^|k=3 , , . . , m} penalty 

i 

n values and one set of assignments for the points in two 

i i 

n 15 dimensions. That is, for each particular there is a 

corresponding two-dimensional penalized cost function problem 
and its solution. Note that the solution of the two- 
dimensional penalized cost function problem identifies the set 
of assignments for the particular values that minimize the 
20 penalized cost function. However, these assignments are not 
likely to be optimum for any higher dimensional problem 
because they were based on an initial arbitrary set of 
values. Therefore, the next step is to determine the optimum 
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assignments for the related three-dimensional penalized cost 
function problem. There exists a two-dimensional hill shaped 
function 0 which is a graph of the minimums of all penalized 
cost functions at various sets of assignments in two 
dimensions. For the three-dimensional problem, the function 
<P can be defined based on the solution to the foregoing two- 
dimensional penalized cost function. By using the current 
values and the {u^|/c > 3} values originally assigned, the 
gradient of the hill -shaped function <P is determined, which 
points toward the peak of the hill. By using the gradient and 
values previously selected for the one point on the hill 
(corresponding to the minimum of the penalized cost function 
^ ) as a starting point, the values can be found for which 
the corresponding problem will result in the peak of the 
function 0. The solution of the corresponding two-dimensional 
problem is the proper values for two of the three sets of 
indices in the three-dimensional problem. These solution 
indices can select a subsection of the cost array which maps 
to a two-dimensional array. The set of indices which minimize 
the two-dimensional assignment problem based on that array 
corresponds to the proper assignment of points in the third 
dimension. The foregoing "recovery" process was known in the 
prior art, but it is modified here to adjust for the sparse 
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matrix characteristic. The next task is to recover the 
solution of the proper values for the four-dimensional 
problem. The foregoing hill climbing process will not work 
again because the foregoing hill climbing process when 
required to locate the four-dimensional values for the peak 
of ^ requires the exact definition of the function <P^ (as was 
available in the case of ^) or an always less than 
approximation of whereas the iteration can result in a 

greater than approximation of <P^ . According to the present 
invention, another three-dimensional function W is defined 
which is a "less than approximation" the three-dimensional 
function 0 and which can be defined based on the solution to 
the two-dimensional penalized cost function and the previously- 
assigned and determined values. Next, the gradient of the 
function W is found and hill climbing technique used to 
determine the values at the peak. Each selected set 
results in a new three-dimensional problem and requires the 
two-dimensional hill climbing based upon new two-dimensional 
problems. At the peak of the three-dimensional function W, 
the solution values are a subset of the values required for 
the four-dimensional solution. Recovery processing extends 
the subset to a complete solution. This process is repeated 
iteratively until the values that result in a corresponding 
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solution at the peak of the highest order functions ¥ and <P 
are found. The final recovery process then results in the 
solution of /c-dimensional problem. 

Figure 5 illustrates process 300 in more detail. 



1.2. Problem Formulation 
In problem formulation step 310, all data structures for 
the subsequent steps are allocated and linked by pointers as 
required for execution efficiency. The incoming problem is 

10 partitioned as described in the subsequent section. This 
partitioning has the effect of dividing the incoming problem 
into a set of independent problems and thus reducing the total 
workload. The partitioning process depends only on the actual 
cost matrix so the partitioning can and is performed for all 

15 levels of the relaxation process. 



Step 32 0 begins the Lagrangian Relaxation process for 
reducing the M-dimensional problem by selecting all Lagrangian 
20 Coefficient penalty values initially equal to zero. The 
Lagrangian Coefficients associated with the constraint set 



dimensional problem in step 322 to a (M-1) -dimensional problem 



5 



1.2.1. Relaxation and Recoverv 



are an (i\7^+l) -element vector. 



The reduction of this M- 
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uses the two step process described above. First, a penalty- 
based on the value of the respective coefficient is added to 
the cost function when a point is assigned to more than one 
track and then the resultant cost function is minimized. 
5 However, during the first iteration, the penalty is zero 
because all values are initially set to zero. Second, the 
requirement that no point from any image can be assigned to 
more than one track is relaxed for one of the images. In the 
extreme this would allow a point from the relaxed image to be 



would not minimize the cost. The effect of the two steps in 
combination is to remove a hard constraint while adding the 



penalty to the cost function so that it operates like a soft 
15 constraint. For step 322 this two-step process results in the 




10 associate with every track. 



previous penalty would probably mean that such an association 



However, the effect of the 



following penalized cost function problem with k^M-. 
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(a) 



0},(u^) = Minimize: (py,(u^,z^) 



i,=0 1,=0 



(b) Subject to: 



(c) 



(1.5) 



(d) 



for 2-,- = 1, , . .N^ 
and j=2, . . .k-1, 

^^...ij, ^ ^0/1^ for all ia...I;c. 



Because the constraint on single assignment of elements from 
the last image has been eliminated, an {M-1) -dimensional 
problem can be developed by eliminating some of the possible 
assignments. As shown in Equations (1.6), this is done by 
selecting the smallest cost element from each of the M^^ axis 
vectors of the cost matrix. Reduction in this manner yields 
a new, lower-order penalized cost function defined by 
Equations (1.6) which has the same minimum cost as does the 
objective function defined by (1.5) (a) above. 

The cost vectors are selected as follows. Define a new 



cost 



k-l 



array c^^...!^.^ by 
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Co. 



'l...i, = Mi-riimum |ci;...i^- ui; | i, = 0, 1, . . . (1.6) 
for {i^, . . . , * {0, . . . ,Q) , 

?. 0 = rnin |o , Co''. . .0 - u J. 
The resulting (M-1) -dimensional problem is (where k.=M) : 
(P^Cu'^; = Minimize: v„.:, ( z"'^ ) = • • • zl c*'.^. .i^_^z*'.^. 



N, N, 



Subject to: £ ■ ■ ■ ll ^i^!'. .i^"^ ' = 1/ • • .^i/ (1-V) 



^2=0 ^;t=0 



5-^ • • • • ■ • 5j ^i, ...i^ 



for = 1, . . .Nj 
and j=2, . . .k-1, 



^i"...! ^{0,1} for all i^...!;,.,, 



Assignment Problem (1.4) and Problem Formulation (1.7) differ 
only in the dimension M vs. M-1, respectively. An optimal 
solution to (1.7) is also an optimal solution to equation 
(1.5). This relationship is the basis for an iterative 
sequence of reductions indicated by steps 320-332 through 330- 
332 and 200-204 in which the penalized cost function problem 
is reduced to a two-dimensional problem. As these formula 
will be used to describe the processing at all levels, the 
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lowercase k is used except where specific reference to the top 
level is needed. In step 206, the two-dimensional penalized 
cost function is solved directly by the prior art Reverse 
Auction technique. Each execution of 206 produces two 
5 results, the set of values that minimize the problem and the 
cost that results when these z values are substituted into 
the objective function (1.7) (a) . 

In step 208, according to the prior art, solution z^ 
values are substituted into the two-dimensional derivative of 

10 the surface The result indicates how the value of 

should be adjusted so as to perform the hill climbing 
function. As was previously described the objective is to 
produce a sequence of values which ends when the Up value is 
in the domain of the peak of the surface <P. The section 

15 "Determining Effective Gradient" describes how new values are 
computed and how it is determined that the current points to 
the peak of When no further adjustment is required the 

flow moves to step 214 which will attempt to recover the 
three-dimensional solution as previously described. When 

2 0 further adjustment is required then the flow progresses to 
step 212 and the new values of are computed. At the two- 
dimensional level the method of the prior art could be used 
for the hill climbing procedure. However, it is not practical 
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to use this prior art hill climbing technique to determine the 
updated Lagrangian Coefficients or the Max on the next (or 
any) higher order surface 0 because the next (or any) higher 
dimensional function 0 cannot be defined based on known 
information. 

Instead, the present invention defines a new function 
based on known information which is useful for hill climbing 
from the third to the fourth and above dimensions, i.e., 
values which result in z values that are closer to the proper 
z values for the highest Jc-dimension . This hill climbing 
process (which is different than that used in the prior art of 
Figure 4 for recovering only the three-dimensional solution) 
is used iteratively at all lower dimensions k of the M- 
dimensional problem (including the three-dimensional level 
where it replaces prior art) even when k is much larger than 
three. Figure 6 helps to explain this new hill climbing 
technique and illustrates the original k-dimensional cost 
function of Problem Formulation (1.4) . However, the actual 
k-dimensional cost surface V;, defined by (1,4) (a) comprises 
scalar values at each point described by ic-dimensional vectors 
and as such can not be drawn. Nevertheless, for illustration 
purposes only, Figure 6 ignores the reality of ordering 
vectors and illustrates a concave function V;c(2^) to represent 
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Equation (1.4) . The surface Vj^ is illustrated as being smooth 
to simplify the explanation although actually it can be 
imagined to be terraced. The goal of the assignment problem 
is to find the values of ~z ^; these values minimize the k- 
5 dimensional cost function Vj^. 

For purposes of explanation, assume that in Figure 6, ^=4 
(the procedure is used for all k ^3) . This problem is reduced 
by two iterations of Lagrangian Relaxation to a two- 

(k-l). 

dimensional penalized cost function 0j . This cost 
10 function, and all other cost functions described below, are 
also non-smooth and continuous but are illustrated in Figure 
6 as smooth for explanation purposes. Solving the (p^ 
problem results in one set of assignments and the value of 

^ (k-l)j^ (k-1). 

^(jc^i). at the point . A series of functions <pj , . . . , (p^ 
15 each generated from a different are shown. The particular 
series illustrates the process of locating the peak ot (P^^^.j^^. 

(k-l). (k-l)j^ 

The two-dimensional penalized cost functions ,...,02 
can be solved directly. Each such solution provides the 
information required to calculate the next value. Each 
20 iteration of the hill climbing improves the selection of 

values, i.e., yields 0^ problem whose solution is closer to 
those at the solution of the 0^. problem. The result of solving 

(k-Dj^ 

01 is values that are on both ^(k-i)x ^k- Figure 6 
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illustrates the surface which comprises the minimums of all 
7c-dimensional penalized cost function surfaces (p^^, i.e., if the 
Problem Formulations (1.5) and (1.7) were solved at all 
possible values of the function ^Pj^iuj^) would result. The 
5 surface is always less than the surface Vj^ except at the 
peak as described in Equation (1.8) and its maximum occurs 
where the surface Vj^ is minimum. Because the surface 
represents the minimum of the surface , any point on the 
surface can mapped to the surface . The function 
10 provides a lower bound on the minimization problem described 
by (1.4) (a) . Let "z ^ be the unknown solution to Problem 
Formulation (1.4) and note that: 



15 Consequently, the values at the peak of (i.e., the 

greatest lower bound on the cost of the relaxed problem) , can 
be substituted into the ^c-dimensional penalized cost function 



present invention attempts to find the maximum on the surface 
2 0 0},. However, it is not possible to use the prior art hill 
climbing to hill climb to the peak of (Pj, because the definition 
of 0^ requires exact knowledge of lower order functions 0. As 
the solution of ^ is not the exact solution, in that higher 



(1.8) 



to determine the proper assignments. 



Consequently, the 
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order values of u are not yet optimized, its value can be 
larger than the true peak of 0^. As such it is not a lower 
bound on v^^ and it can not be used to recover the required 
solution. 

Instead, the present invention defines all auxiliary 
functions which are based only on the solution to the two- 
dimensional penalized cost function problem, lower order 
values and values determined previously by the reduction 
process. The function W^^ is a less than approximation of 0^^, 
and its gradient is used for hill climbing to its peak. The 
values at the peak of the function are then substituted 
into Problem Formulation (1.5) to determine the proper 
assignments. To define the function W^, the present invention 
explicitly makes the function ^;,(u^) a function of all higher 
order sets of Lagrangian Coefficients with the expanded 
notation: 0i^[u^ ; u^""^, . . . , u^) . Then, a new set of functions W, 
is defined recursively, using the 0)^' s domain: 

W3(u3) = V2 +5^ =03(u^•u^ . . , (1.9) 

i3 = 0 ' 

where V2 is the solution value for the most recent two- 
dimensional penalized cost function problem. For k >3 
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if known, 



(1.10) 



othert/ise. 



'43 



From the definitions of and V;, (Problem Formulation (1,7) 
compared with Problem Formulation (1.5)) : 



(Dj^ ( u u . . . , u ^) = v^., ( z ) + 2J u/^. 



it follows that : 



and with that Equation (1.8) is extended to: 

w^(u^ . . .,u^-^•u^)^<^>^(u^•u^^^ . . .,u^)^ v^(u^)^ v^(u^) (i.id 

This relationship means that either 0^ or W^. may be used in 
hill climbing to update the Lagrangian Coefficients u. <t>k 
the preferred choice, however it is only available when the 
solution to Problem Formulation (1.5) is a feasible solution 
to Problem Formulation (1.4) (as in hill climbing from the 
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second to third dimension which is why prior art worked) . For 
simplicity in implementation, the function W^^ is defined so 
that it equals (Pj^ when either function could be used. It is 
therefore always used, even for hill climbing from the second 
5 dimension. The use of the function W which is based on 
previously assigned or determined u values and not higher 
order u values which are not yet determined, is an important 
feature of the present invention . 

10 1.2.2. Determining Effective Gradient 

After the function 5^ is defined, the next steps of hill 
climbing/peak detection are to determine the gradient of the 
function W^, determine an increasing portion of the gradient 
and then move up the surface Wy, in the direction of this 

15 increasing portion of the gradient. As shown in Figure 6, any 
upward step on the surface W^^, for example, to the minimum of 
the (p^^ will yield a new set of values (to the "left") that 
is closer to the ideal set of values which correspond to the 
minimum of the function v^. While it is possible to 

20 iteratively step up this surface with steps of fixed size 
and then determine if the peak has been reached, the present 
invention optimizes this process by determining the single 
step size from the starting point at the minimum of that 
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will jump to the peak and then calculate the values at the 



peak. Once the values u at the peak of Wj, are determined, then 



the values can be substituted into Problem Formulation (1.5) 



to determine the proper assignment. (However, in a more 



5 realistic example, where k is much greater than three, then 



the values at the peak of the function ¥ along with the 



lower-order values xsl^ and those assigned and yielded by the 



reduction steps are used to define the next higher level 



t function ¥. This definition of a higher order function ¥ and 

f] 10 hill climbing process are repeated iteratively until ¥^, the 



peak of ¥i^f and the values at the peak of ¥j^ are identified.) 



The following is a description of how to determine the 



gradient of each surface ¥ and how to determine the single 



step size to jump to the peak from the starting point on each 



15 surface ¥. 



As noted above, each surface ¥ ±s non-smooth. Therefore, 



if a gradient was taken at a single point, the gradient may 



not point toward the peak. Therefore, several gradients (a 



"bundle") are determined at several points on the surface if^, 



20 in the region of the starting point (i.e., minimum of 0^,^) and 



then averaged. Statistically, the averaged result should 



point toward the peak of the surface ¥^. Wolfe's Conjugate 



Subgradient Algorithm (P. Wolfe. A method of conjugate 
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subgradients for minimizing non-dif f erentiable functions. 
Mathematical Programming Study, 3:145-173, 1975; P. Wolfe. 
Finding the nearest point in a polytope. Mathematical 
Programming Study, 11:128-149, 1976) for minimization was 
5 previously known in another environment to determine a 
gradient of a non-smooth surface using multiple subgradients 
and can be used with modification in the present invention. 
The modification to Wolfe's algorithm uses the information 
generated for }Fj^{u^^, . . .,u^^"^; u^^"^) as the basis for calculating 
10 the subgradients. The definition of a subgradient v of 
^'^.(u^^, . . .,u^^) is any member of the subdif f erential set defined 
as : 



15 ^v^(u'-u^) V u'G R^J^^^}, 

where is the transpose of v. 

Next, a subgradient vector is determined from this 

function. If is the solution of Problem Formulation (1.5), 

then differentiating !f^;,(u^, . . . ,u^"^; u*^) with respect to Uij, and 
20 evaluating the result with respect to the current selection 

matrix yields a subgradient vector: 
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g= (0, (1 



-t...t 



= 1. 



The iterative nature of the solution process at each dimension 
yields a set of such subgradients . Except for a situation 
described below, where the resultant averaged gradient does 
not point toward the peak, the most recent set of such 
subgradients are saved and used as the "bundle" for the peak 
finding process for this dimension. For example, at the 
level k there is a bundle of subgradients of the surface 
near the minimum of the surface (p^. determined as a result of 
solving Problem Formulation (1.4) at all lower levels. This 
bundle can be averaged to approximate the gradient . 
Alternately, the previous bundle can be discarded so as to use 
the new value to initiate a new bundle. This choice provides 
a way to adjust the process to differing classes of problems, 
i.e., when data is being derived from two sensors and the 
relaxation proceeds from data derived from one sensor to the 
other then the prior relaxation data for the first sensor 
could be detrimental to performance on the second sensors 
data . 
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1.2.3. Step Size and Termination Criterion 
After the average gradient of the surface ¥^ is 
determined, the next step is to determine a single step that 
will jump to the peak of the surface W^. The basic strategy 
is first to specify an arbitrary step size, and then calculate 
the value of at this step size in the direction of the 
gradient. If the value of is larger than the previous one, 
this probably means that the step has not yet reached the 
peak. Consequently, the step size is doubled and a new value 
of fjj. is determined and compared to the previous one. This 
process is repeated until the new value of W^^ is less than the 
previous one. At that time, the step has gone too far and 
crossed over the peak. Consequently, the last doubling is 
rolled-back and that step size is used for the iteration. If 
the initial estimated is less than previous value then the 
step size is decreased by 50%. If this still results in a 
smaller value of iPj,, then the last step is rolled back and the 
previous step size is decreased by 25%. The following is a 
more detailed description of this process. 

With a suitable bundle of subgradients determined as just 
described, Wolfe's algorithm can be used to determine the 
effective subgradient d and the upgraded value Uj^^. From the 
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previous iteration, or from an initial condition, there exists 
a step length value t. The value, 

= Uj + td 

is calculated as an estimate of Uj+i. To determine if the 
5 current step size is valid we evaluate ^^^(u^, . . . /U^'^; u^) - If 
the result represents an improvement then we double the step 
size. Otherwise we halve the step size. In either case a new 
is calculated. The doubling or halving continues until the 
step becomes too large to improve the result, or until it 
10 becomes small enough to not degrade the result. The resulting 
suitable step size is saved with d as part of the subgradient 
bundle. The last acceptable u^. is assigned to u^^^ . 

Three distinct criteria are used to determine when Uj is 
close enough to u^: 
15 1. The Wolfe's algorithm criterion of d=0 given that the 
test has been repeated with the bundle containing only 
the most recent subgradient. 

2. The difference between the lower bound ^j.Cu'^) and the 
upper bound v^Cz^", u^) being less than a preset relative 

2 0 threshold. (Six percent was found to be an effective 

threshold for radar tracking problems.) 

3. An iteration count being exceeded. 
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The use of limits on the iteration are particularly 
effective for iterations at the level 3 through (n-1) as these 
iterations will be repeated so as to resolve higher order 
coefficient sets. With limited iterations the process is in 
5 general robust enough to improve the estimate of upper order 
Lagrangian Coefficients. By limiting the iteration counts 
then the total processing time for the algorithm becomes 
deterministic. That characteristic means the process can be 
effective for real time problems such as radar, where the 
10 results of the last scan of the environment must be processed 
prior to the next scan being received. 

1.2.4. Solution Recovery 
The following process determines if the values at what 

15 is believed to be the peak of the function adequately 
approximate the values at the peak of the corresponding 
function 0^^. This is done by determining if the corresponding 
penalized cost function is minimized. Thus, the values at 
the peak of the function are first substituted into Problem 

2 0 Formulation (1.5) to determine a set of z assignments for k~l 
points. During the foregoing reduction process. Problem 
Formulation (1.7) yielded a tentative list of k selected z 
points that can be described by their indices as: 
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. . • / where Nq is the number of cost elements selected 
into the solution. One possible solution of Problem 
Formulation (1.4) is the solution of Problem Formulation (1.5) 
which is described as ii . - - i^-i^i^. . J j with as it 

was defined in Problem Formulation (1.6). If this solution 
satisfies the constraint set, then it is the optimal 

solution. 

However, if the solution for Problem Formulation (1.5) is 
not feasible (decision 355) , then the following adjustment 
process determines if a solution exists which satisfies the k^^ 
constraint while retaining the assignments made in solving 
Problem Formulation (1.7) . To do this, a two-dimensional cost 
matrix is defined based upon all observations from the k^^ set 
which could be used to extend the relaxed solution. 

for 1 = 0, . . . ,Nj^ 
^ji^^^i Vi'-^ and j = 0, . . . ,N^. (1.12) 

If the resulting two-dimensional assignment problem. 
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Minimize 
Subject to: 



5 

has a feasible solution, then the indices of that solution map 
to the solution of Problem Formulation (1.4) for the k- 
dimensional problem. The first index in each resultant 
solution entry is the pointer back to an element of the 

10 1^ if'. . . i^_^in^^_ ..^ J I list. That element supplies the first k-1 
indices of the solution. The second index of the solution to 
the recovery problem is the k^^ index of the solution. 
Together these indexes specify the values of that solve 
Problem Formulation (1.4) at the k^^ level. 

15 If Problem Formulation (1.13) does not have a feasible 

solution, then the value of Ui which was thought to represent 
Up is not representative of the actual peak and further 
iteration at the k^ level is required. This decision 
represent the branch path from step 214 and equivalent steps. 

20 



j = 0 1 = 0 ^ 

i; w =1 

1=0 ■' 



j=l, . . . ,No, 



(1.13) 



1=0 -' 

w 6(0,1} j = o, 



.,No, 1 = 0, 
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1.3. Partitioning 



A partitioning process is used to divide the cost matrix 
that results from the Scoring step 154 into as many 
independent problems as possible prior to beginning the 
5 relaxation solution. The partitioning process is included 



partitioning is a set of problems to be solved, i.e., there 
will be Pi problems that consist of a single hypothesis, Pz 
problems that consist of two hypothesis, etc. Each such 

10 problem is a group in that one or more observations or tracks 
are shared between members of the group. 

In partitioning to groups no consideration is given to 
the actual cost values. The analysis depends strictly on the 
basis of two or more cost elements sharing the same specific 

15 axis of the cost matrix. In a two-dimensional case two cost 
elements must be in the same group if they share a row or a 
column. If the two elements are in the same row, then each 
other cost element that is also in the same row, as well as 
any cost elements that are in columns occupied by members of 

2 0 the row must be included in the group. The process continues 
recursively. In literature it is referred to as "Constructing 
a Spanning Forest." The ic-dimensional case is analogous to 
the two-dimensional case. The specific method we have 



with the Problem Formulation Step 310. 



The result of 
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incorporated is a depth first search, presented by Aho, 
Hopecraft, and Ullman (A.V. Aho, J.E. Hopcroft, and J.D, 
Ullman. Design and Analysis of Computer Algorithms , Addi son- 
Wesley, MA, 1974) . 

The result of partitioning at level M is the set of 
problems described as {P^j| i=l , . . . ,Pj- and j =1 , . . . , i\7} , where N is 
the total number of hypothesis. The problems labeled 
{Pii\i=l , ' ' ' , Pi} are the cases where there is only one choice 
for the next observation at each scan and that observation 
could be used for no other track, i.e., it is a single 
isolated track. The hypothesis must be included in the 
solution set and no further processing is required. 

As hypothesis are constructed the first element is used 
to refer to the track ID. Any problem in the partitioned set 
which does not have shared costs in the first scan represent 
a case where a track could be extended in several ways but 
none of the ways share an observation with any other track. 
The required solution hypothesis for this case is the 
particular hypothesis with the maximum likelihood. For this 
case all likelihoods are determined as was described in 
Scoring and the maximum is selected. 

In addition to partitioning at the level M, partitioning 
is applied to each subsequent level M-l,...,2. For each 
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problem that was not eliminated by the prior selection, 
partitioning is repeated, ignoring observations that are 
shared in the last set of observations. Partitioning 
recognizes that relaxation will eliminate the last constraint 
5 set and thus partitioning is feasible for the lower level 
problems that will result from relaxation. This process is 
repeated for all levels down to k^2 , The full set of 
partitionings can be performed in the Problem Formulation Step 
310, prior to initiating the actual relaxation steps. The 

10 actual two-dimensional solver used in step 206 includes an 
equivalent process so no advantage would be gained by- 
partitioning at the k=2 level. 

There are two possible solution methods for the remaining 
problems. "Branch and Bound" as was previously described, or 

15 the relaxation method that this invention describes. If any 
of the partitions have 5-10 possible tracks and less than 50 
to 2 0 hypotheses, then the prior art "Branch and Bound" 
algorithm generally executes faster than does the relaxation 
due to its reduced level of startup overhead. The "Branch and 

2 0 Bound" algorithm is executed against all remaining M level 
problems that satisfy the size constraint. For the remainder 
the Relaxation algorithm is used. The scheduling done in 
Problem Formulation allows each Problem Formulation (1.5) cost 



CSUR.01USRI 

matrix resulting from the first step of relaxation to be 
partitioned. The resulting partitions can be solved by any of 
the four methods: isolated track direct inclusion, isolated 
track tree evaluation, small group "Branch and Bound" or an 
additional stage of relaxation as has been fully described. 

The partitioning after each level of Lagrangian 
Relaxation is effective because when a problem is relaxed, the 
constraint that requires each point to be assigned to only one 
track is eliminated (for one image at a time) . Therefore, two 
tracks previously linked by contending for one point will be 
unlinked by the relaxation which permits the point to be 
assigned to both tracks. The fruitfulness of partitioning 
increases for lower and lower levels of relaxation. 

The following is a more detailed description of the 
partitioning method. Its application at all but the level M 
depends upon the relaxation process described in this 
invention. The recursive partitioning is therefore a distinct 
part of this invention. The advantage of this method is 
greatly enhanced by the sparse cost matrix resulting from 
tracking problems. However the sparse nature of the problem 
requires special storage and search techniques. 

A hypothesis is the combination of the cost element the 
selection variable and all observations that made up the 
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potential track extension. It can be written as, 

K^[c^,z^, { 0,^=0;,. I ;c=l, . . . E {1, . . . / N;c}}}/ i-e. , cost, 

selection variable and observations in the hypothesis. Here 
n €{l, . . . ,]<!], where N is the total number of hypotheses in the 
problem. While the cost and assignment matrices were 
previously referenced, these matrices are not effective 
storage mechanisms for tracking applications. Instead the 
list of all hypothesis and sets of lists for each dimension 
that reference the hypothesis set are stored. The 
hypothetical set in list form is: 

For each dimension Jc=l,...,M there exists a set of lists, 
with each list element being a pointer to a particular 
hypothesis : 

i=Wj^ ,k=M 

where N^. is a number of hypothesis containing the i^^ 
observation from scan k. This structure is a multiply linked 
list in that any observation is associated with a set of 
pointer to all hypothesis it participates in, and any 
hypothesis has a set of pointers to all observations that 
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formed it. (These pointers can be implemented as true 
pointers or indices depending upon the particular programming 
language utilized.) 

Given this description of the matrix storage technique 
5 then the partitioning technique is as follows: Mark the first 
list L^^ and follow out all pointers in that list to the 
indicated hypothesis hp^^ for i=l,.,.,N^ . Mark all located 
hypothesis, and for each follow pointers back the particular 
for k=\, , . . ,M. Those L' s if not previously marked get 

10 marked and also followed out to hypothesis elements and back 
to L*s. When all such L's or h's being located are marked, 
then an isolated sub-problem has been identified. All marked 
elements can be removed from the lists and stored as a unique 
problem to be solved. The partitioning problem then continues 

15 by again starting at the first residual set L. When none 
remain, the original problem is completely partitioned. 

Isolated problems can result from one source track having 
multiple possible extensions or from a set of source tracks 
contending for some observations. Because one of the indices 

20 of k (in our implementation it is k=l) indicates the source 
track then it is possible to categorize the two problem types 
by observing if the isolated problem includes more than one 
L-list from the index level associated with tracks. 
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1.4. Subsequent Track Prediction 
As noted above, after the points are assigned to the 
respective tracks, some action is taken such as controlling 
take-offs and landings to avoid collision, advising airborne 
5 aircraft to change course, warning airborne aircraft of an 
adjacent aircraft in a commercial environment, or aiming and 
firing an anti-aircraft (or anti-missile) missile, rocket or 
projectile, or taking evasive action in a military 
environment. Also, the tracking can be used to position a 

10 robot to work on an object. For some or all of these 
applications, usually the tracks which have just been 
identified are extended or extrapolated to predict a 
subsequent position of the aircraft or other object. The 
extrapolation can be done in a variety of prior art ways; for 

15 example, straight line extensions of the most likely track 
extension hypothesis, parametric quadratic extension of the x 
versus time, y versus time, etc. functions that are the result 
of the filtering process described earlier, least square path 
fitting or Kalman filtering of just the selected hypothesis. 

20 The process of gating, filtering, and gate generation as 

they were previously described require that a curve be fit 
through observations such that fit of observations to the 
likely path can be scored and further so that the hypothetical 
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target future location can be calculated for the next stage of 
gating. In the implementation existing, quadratic functions 
have been fit through the measurements that occur within the 
window. A set of quadratics is used, one per sensor 
5 measurement. To calculate intercept locations these 

quadratics can be converted directly into path functions. 
Intercept times are then calculated by prior methods based 
upon simultaneous solution of path equations. 

The use of the fitted quadratics is not as precise as 

10 more conventional filters like the Kalman Filter. They are 
however much faster and sufficiently accurate for the relative 
scores required for the Assignment problem. When better 
location prediction is required then the assignment problem is 
executed to select the solution hypothesis and based upon the 

15 observations in those hypothesis the more extensive Kalman 
filter is executed. The result is tremendous computation 
savings when compared with the Kalman Filter being run on all 
hypothetical tracks . 

Based on the foregoing, apparatus and methods have been 

20 disclosed for tracking objects. However, numerous 

modifications and substitutions can be made without deviating 
from the scope of the present invention. For example, the 
foregoing functions W can also be defined as recursive 
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approximation problem in which several values of higher order 
u'^ values are used to eliminate the higher than approximation 
characteristic of the function 0^^, The hill climbing of the 
function W can be implemented by a high order hill climbing 
5 using the enhanced function <P^. Although the result would not 
be as efficient it seems likely that the method would 
converge. Therefore, the invention has been disclosed by way 
of example and not limitation, and reference should be made to 

r1 the following claims to determine the scope of the present 

'4 J 

i;0 10 invention. 

fn 

II. AN ALTERNATIVE EMBODIMENT OF THE MULTI -DIMENSIONAL 
ASSIGNMENT SOLVING PROCESS 

Q The following description provides an alternative second 

U I 

ry 15 embodiment of the multi -dimensional assignment solving process 

U1 

Q of section I.l. However, in order to clearly describe the 

present alternative embodiment, further discussion is first 
presented regarding the data assignment problem of 
partitioning measurements into tracks and false alarms and, in 
2 0 particular, the representation of multi-dimensional assignment 
problems . 
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II. 1. Formulation of the Assignment Problem 
The goal of this section is to explain the formulation of 
the data association problems, and more particularly multi- 
dimensional assignment problems, that govern large classes of 
problems in centralized or hybrid centralized-sensor level 
multisensor/multitarget tracking. The presentation is brief; 
technical details are presented for both track initiation and 
maintenance in (A.B. Poore . Multidimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:21-51, 1994) 
for non-maneuvering targets and in (A.B. Poore. 
Multidimensional assignments and multitarget tracking: 
Partitioning data sets. In P. Hansen, I.J. Cox, and B. 
Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, volume 19, pages 169-198, 
Providence, R.I., 1995. American Mathematical Society) for 
maneuvering targets. Note that the present formulation can 
also be modified to include target features (e.g., size and 
type) into the scoring process 154 . 

The data assignment problems for multisensor and 
multitarget tracking considered here are generally posed as 
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that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



Maximize 



P(r=Y°|Z^) 



(2.1) 



where represents M data sets, y is a partition of indices 
of the data (and thus induces a partition of the data) , F* is 
the finite collection of all such partitions, T is a discrete 
random element defined on T*, y° is a reference partition, and 
10 P{T=y\z^) is the posterior probability of a partition y being 
true given the data Z^. The term partition is defined below. 

Consider M observation sets Z (A:) (A:=l , . . . , i\7) each of 
observations ] , and let denote the cumulative data set 
defined by 

Z{k) and Z'" = {z(l) , . . . , Z(M)}, (2.2; 



15 respectively. In multisensor data fusion and multitarget 
tracking the data sets Z{k) may represent different classes of 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must be 
partitioned into tracks and false alarms. In the formulation 

2 0 of track extensions a moving window over time of observations 
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sets is used. The observation sets will be measurements which 
are: (a) assigned to existing tracks, (b) designated as false 
measurements , or (c) used for initiating tracks . However, 
note that alternative data objects instead of observation sets 
5 may also be fused such as in sensor level tracking wherein 
each sensor forms tracks from its own measurements and then 
the tracks from the sensors are fused in a central location. 
Note that, as one skilled in the art will appreciate, both 
embodiments of the present invention may be used for this type 

10 of data fusion. 

The data assignment problem considered presently is 
represented as a case of set partitioning defined in the 
following way. First, for notational convenience in 

representing tracks, a zero index is added to each of the 

15 index sets in (2.2), a gap filler is added to each of the 
data sets Z{k) in (2.2), and a "track of data" is defined as 



partition of the data will refer to a collection of tracks of 
data or track extensions wherein each observation occurs 
2 0 exactly once in one of the tracks of data and such that all 
data is used up. Note that the occurrence of the gap filler 
is unrestricted. The gap filler serves several purposes in 
the representation of missing data, false observations, 




where i^^ can now assume the values of 0 . A 
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initiating tracks, and terminating tracks. Note that a 
"reference partition" may be defined which is a partition in 
which all observations are declared to be false. 

Next under appropriate independence assumptions it can be 
shown that : 



-^llrdZ^.L. n L. . , (2.3) 



-1- • - -W 



where Lf^...!^ is a likelihood ratio containing probabilities for 
10 detection, maneuvers, and termination as well as probability 
density functions for measurement errors, track initiation and 
termination. Then if Ci^...ij^= -In (L^^^ .i^) , 



15 



-In 



P(r = Y°|2^ 



(1^, . . .,i^)eY 



(2.4) 



Using (2.3) and the zero-one variable 2^^,..^^=!, if 

E Y/ ^rid 0, otherwise, one can then write the 
problem (2,1) as the following M-dimensional assignment 
problem (as also presented in section Problem Formulation 
20 (1.4) with k=M) : 
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(a) Minimize 



(b) SuJbject to: 



1,-0 vo -i---^" 



(c) 



E ... S S (2-^^ 

for i;, = 1, . . . and k^2, . . . ,M-1, 
5 (d) E • • • E i =1/ = ^/ • • • '^M/ 

(e) ^ ^O'l^ for all i,, 

where c^..-^ is arbitrarily defined to be zero. Here, each 
group of sums in the constraints (2.5) (b)-(2.5) (e) represents 

10 the fact that each non-gap filler observation occurs exactly 
once in a "track of data." One can modify this formulation to 
include multi-assignments of one, some, or all the actual 
observations. The assignment problem (2.5) is changed 
accordingly. For example, if is to be assigned no more 

15 than, exactly, or no less than n)^^ times, then the " = 1" is 
changed to =, ^ ^i;^' " respectively. 

Expressions for the likelihood ratios L\^,_i^ appearing in 
the costs c. . =-ln(L. . ) are well-known. In particular, 
discussions of these ratios may be found in (A.B. Poore. 

20 Multidimensional assignment formulation of data association 
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problems arising from multitarget tracking and multisensor 
data fusion, Computational Optimization and Applications , 
3:27-57, 1994; A.B. Poore . Multidimensional assignments and 
multitarget tracking: Partitioning data sets. In P. Hansen, 
5 I.J. Cox, B. Julesz, editors, DIMACS Series in Discrete 
Mathematics and Theoretical Computer Science, volume 19, pages 
169-198, Providence, R.I., 1995. American Mathematical 
Society) . Furthermore, the likelihood ratios are easily 
modified to include target features and to account for 

10 different sensor types. Also note that in track initiation, 
the M observation sets provide observations from M sensors, 
possibly all the same. Additionally note that for track 
maintenance, a sliding window of M observation sets and one 
data set containing established tracks may be used. In this 

15 latter case, the formulation is the same as above except that 
the dimension of the assignment problem is now M+1. 



II. 2. Overview of the New Lagrancrian 
Relaxation Algorithms 

2 0 Having formulated an M-dimensional assignment problem 

(2.5), we now turn to a description of the Lagrangian 
relaxation algorithm. Subsection II. 2.1 below presents many 
of the relaxation properties associated with the relaxation of 
an n-dimensional assignment problem to an in-dimensional one 
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via a Lagrangian relaxation of n-m sets of constraints wherein 
m < n < M and preferably in the present invention embodiment 
n-m >1. Although any n~m constraint sets can be relaxed, the 
description here is based on relaxing the last n-m sets of 
constraints and keeping the first m sets. Given either an 
optimal or suboptimal solution of the relaxed problem, a 
technique for recovering a feasible solution of the 
n-dimensional problem is presented in subsection II. 2. 2 below 
and an overview of the Lagrangian relaxation algorithm is 
given in subsection II. 2. 3. 

The following notation will be used throughout the 
remainder of the work. Let M be an integer such that M >3 and 
let nE{3,...,M}. The n-dimensional assignment problem is 



(a) Minimize: v^('zj 




(b) Subject To 




(2.6) 



(c) 




for i,=l 



and k = 2 



f ' • • f 



n-l, 



(d) 




(e) 



P V elo, l] for all i 

1 • ' n ^ * 
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To ensure that a feasible solution of (2.6) always exists, all 
variables with exactly one nonzero index (i.e., variables of 
the form z^. . .Qi^Q. . for i^^^Q) are assumed free to be assigned 
and the corresponding cost coefficients are well-defined. 

5 

II. 2.1. The Laqranqian Relaxed Assignment Problem 
The n-dimensional assignment problem (2.6) has n sets of 
constraints. A (i\/i,+ l) -dimensional multiplier vector associated 
with the k^^ constraint set will be denoted by u^=(Uo,U2, . . . , u^^) ^ 

10 with Uo=0 and A:=l,...,n. The n-dimensional assignment problem 
(2.6) is relaxed to an m-dimensional assignment problem by 
incorporating n-m of the n sets of constraints into the 
objective function (2.6) (a) . Although any constraint set can 
be relaxed, the description of the relaxation procedure for 

15 (2.6) will be based on the relaxation of the last n-m sets of 
constraints. The relaxed problem is 

0^^(u^"S . . . , u") = Minimize ( z ^; u^'S . . . , u") (2.7) 

\^ /I 

* 1 1 < 



•'• 2^ • - • ^i, 



1 
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Subject To 



t. 


..i 






t. 









k 

r . 



for = 1, . . . , and A: = 2, . . . ,m, 
2v..i„^{0^l} ^or all 2^, . . . , i^. 



wherein we have added the multiplier Uo=0 for notational 
convenience. Thus, the multiplier u^eR^''^'^ with Uo=0 is fixed. 

An optimal (or suboptimal) solution of (2.7) can be 
constructed from that of an xn-dimensional assignment problem. 
To show this, define for each (i-^, . . . , i^) an index (j^n+i, . . . ,jn) 
= (Jin+i (ii/ — . / in,) / • • • / jn(ii/ • • • / iin) ) ^nd a new cost function 



= argmin ^ 



W. and k = m-^l, . . . , n} f (2.8) 



:f i =c. • + Lj-f -for (i-,...,i 



ij ^ (0, . . .0) , 
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If ( Jin+i/ • • • / jn) is not unique, choose the smallest such j^^^, 
amongst those {n-m) -tuples with the same j^^j. choose the 
smallest j^^2/ etc., so that ( j^n+i/ • • • / jn) is uniquely defined.) 
Using the cost coefficients defined in this way, the following 
/n-dimensional assignment problem is obtained: 
^m.(u"''/ . • wu") = Minimize (2^;u^^\ . . . , u^) = v^(z^) 

il=0 i„=0 1 » 1 " 

Subject to: • • • X/ ...i =1/ -ii = 1/ • • • (2 . 9) 



i2 = 0 -i™=0 



5^ • • - S ^ • - - 52 i' 

i,=o v,=o i^,,=o i„=o 
for i;^. = 1, . . . ,Ny. and k=2, . . . ,m-l, 



• • • -^2, ... 1 i / i/n"^ / • • • / -^m/ 

r,...i, ^ ^0.1^ for all i,, . . 



As an aside, observe that any feasible solution z"^ of (2.6) 
yields a feasible solution of (2.9) via the construction 



m 

1 tt 



0, otherwise. 
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Thus, the jn-dimensional assignment problem (2.9) has at least 
as many feasible solutions of the constraints as the original 
problem (2.6) . 

The following Fact A. 2 has been shown to be true. It 
5 states that an optimal solution of Problem Formulation (2.7) 
can be computed from that of Problem Formulation (2.9). 
Furthermore, the converse to this fact is provided in Fact 
A. 3. Moreover, if the solution of either of these two 
problems (2.7) or (2.9) is E-optimal (i.e., the objective 
ff\ 10 function associated with the suboptimal solution is within "G" 

: 

4' of the objective function associated with the optimal 

ijl solution) , then so is the other. 

□ Fact A. 2. Let be a feasible solution to problem (2.9) 

and define w" by 



(ij, . . .,ij*{0, . . .,0) 

w2^...i^=0 if . . .i„) " . . .Jn) and (2.10) 

(ii, . . . # (0, . . . ,0) 



^ oi , i =0 i-f cS oi . i + /J >Q' 



Then w" is a feasible solution of the Lagrangian relaxed 
probl em (2.7) and 

86 



CSUR.01USRI 



,m+l 



,u") = 



n+1 



f • • • f 



£ t < 



If, in addition, is optimal for (2.9), then is an optimal 
solution of (2.7) and 

^^(u-^...,u") = ^l(u-^S...,u")- £ 1^ u," . 

Jt=2n+1 ij^=0 



With the exception of one equality being converted to an 
inequality, the following Fact is a converse of Fact A. 2 and 
has also been shown to be true. 

Fact A. 3. Let w" be a feasible solution to problem (2.7) 
and define by 

^i,..in.= zl w?,...i, for (i,, . . .,i,)^(0, 0), 

wS...o=0 if (ii, . . .,iJ-(0, . . . ,0) and 

^...ow...i;, + £ "i", >0 for all ' ' ' fin) i (2.11) 

wS...o=l if (i:i, . . .,in,) = (0/ - . .,0) and 
^...ow...i, + £ "il ^0 for some ( 

;c=m+l 
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Then vf ±s a. feasible solution of the problem (2.9) and 



If, in addition, W is optimal for (2.7) , then is an optimal 
solution of (2.9), 

k=m+l ij^=0 
k 



II. 2. 2, The Recovery Procedure 
The next objective is to explain a recovery procedure or 
method for recovering a solution to the n -dimensional problem 
of (2.6) from a relaxed problem having potentially 
substantially fewer dimensions than (2.6). Note that this 
aspect of the alternative embodiment of the present invention 
is substantially different from the method disclosed in the 
first embodiment of the multidimensional assignment solving 
process of section I.l of this specification. Further, this 
alternative embodiment provides substantial benefits in terms 
of computational efficiency and accuracy over the first 
embodiment, as will be discussed. Thus, given a feasible 
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(optimal or suboptimal) solution u/" of (2.9) (or if Problem 
Formulation (2.7) is constructed via Fact A. 2), an explanation 
is provided here regarding how to generate a feasible solution 
of (2.6) which is close to in a sense to be specified. 
5 We first assume that no variables in (2.6) are preassigned to 
zero (this assumption will be removed shortly) . The 
difficulty with the solution w" is that it need not satisfy the 
last n-m sets of constraints in (2.6) . Note, however, that if 
ijfj is an optimal solution for (2.9) and (constructed as in 

jjl 10 Fact A. 2) satisfies the relaxed constraints, then w" is optimal 
for (2.6) . The recovery procedure described here is designed 
to preserve the 0-1 character of the solution of (2.9) as 
far as possible. That is, if w^^...i^=l and ij^O for at least one 
1=1, . . . ,n3, then the corresponding feasible solution z"" of (2.6) 
15 is constructed so that Zi^...i^=l for some (i^^^, . . . , i^,) . Note 
that by this reasoning, variables of the form 2o. . .oi^n+j. . .i^^ tie 
assigned to a value of 1 in the recovery problem only if 
wS_,o=l. However, variables ^'o,..oi^^i,..in be treated 

differently in the recovery procedure in that they can be 
2 0 assigned 0 or 1 independent of the value wS...o- This increases 
the feasible set of the recovery problem, leading to a 
potentially better solution. 



in 
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Let I (i/, . . . , il) 1*'°^ be an enumeration of indices of w™ (or 
the first m indices of w" constructed as in Fact A. 2) such that 
w"^^ .,. = 1 and (ii, . . ^ (0, . . . ,0) . Set . . . , i°) = (0, . . . , 0) 

for j-0 and define 



i = <^.j .J. . 2 , 12 

for i;c=0/ • - • t^ki ^=-m+l, . . . ,n; j=0, . . . 



Let Y denote the solution of the (n-jn+l) -dimensional 
10 assignment problem: 



Minimize 2I zl * — zl ^ir\\ ,i Yn i 



N.., N, 



Subject to ^ • • • 5^ y-ii i =1' j=l,...,No, (2.13) 



2^ 2^ ... 2^ 2^ ...l^y^i ,...i=i/ 

j = 0 i„,i = 0 Vi = 0 l,,i = 0 i„=0 -1 " 



for i^=l, . . . and k=in+l, . . . ,n-l, 

15 22 Z2 " ' Zl Yii i =1' in=I/ • . • /-Nm 



The recovered feasible solution z" of (2.6) corresponding to 
20 the multiplier set {u", ...,u"} is then defined by 
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n 



1, if [i^, . . . , ij = [il, . . . , i^) for some 

j^O, . . . ,N. and Y., ...i=l, (2.14) 

0, otherwise. 



This recovery procedure is valid as long as all cost 
coefficients are defined and all zero-one variables in 
are free to be assigned. Note that modifications are 
necessary for sparse problems. If the cost coefficient 
5 c^j .is well defined and the zero-one variable 

.j, .is free to be assigned to zero or one, then define 
c--''"^ i =c^ J as in (2.12) with zj^""'^ ^ being free to 

be assigned. Otherwise, ^jV^^."^. .i preassigned to zero. To 

J zn+l • * * n 

ensure that a feasible solution exists, we now only need 
10 ensure that the variables Zjo^^^o are free for j=0 , 1 , . . . , . 
(Recall that all variables of the form z^.^^^^.o are free (to be 
assigned) and all coefficients of the form c?...!^. ..o well 
defined for k=l,.,,,n,) This is accomplished as follows. If 
the cost coefficient c". . . is well defined and 

15 z". is free, then define c^n"^''\-c^i ^ with z!ln^^\ 

being free. Otherwise, since all variables of the form 
^o...ii,...o known feasible and have well-defined costs, put 



n-m+l _ X ^ n 
CjO...O- ^O...0i;/O...0 
A:=l, ij; *o 
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II. 3. The Mult i -Dimensional LaQranqian Relaxation 

Algorithm for the n-Dimensional Assignment 
Problem 

Starting with the M-dimensional assignment problem (2.6) , 
5 i.e., n^M, the algorithm described below is recursive in that 
the M-dimensional assignment problem is relaxed to an m- 
dimensional problem by incorporating {n-m) sets of constrains 
into the objective function using Lagrangian relaxation of 
this set. This problem is maximized with respect to the 

10 Lagrange multipliers, and a good suboptimal solution to the 
original problem is recovered using an (n-in+1) -dimensional 
assignment problem. Each of these two (the m-dimensional and 
the (n--m+l) -dimensional assignment problems) can be solved in 
a similar manner. 

15 More precisely, reference is made to Figure 8 which 

presents a flowchart of a procedure embodying the multi- 
dimensional relaxation algorithm, referred to immediately 
above. This procedure, denoted MULTI_DIM_REIiAX in Figure 8, 
has three formal parameters, n, PROB_FORMUIiATION and 

20 ASSIGNMT_SOLUTION, which are used to transfer information 
between recursive instantiations of the procedure. In 
particular, the parameter, n, is an input parameter supplying 
the dimension for the mult i -dimensional assignment problem (as 
in (2. 6)) which is to be solved (i.e., to obtain an optimal or 
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near -optimal solution) . The parameter, PROB_FORMULATION, is 
also an input parameter supplying the data structure (s) used 
to represent the n-dimensional assignment problem to be 
solved. Note that PROB_FORMULATION at least provides access 
5 to the cost matrix, [c"] , and the observation sets whose 
observations are to be assigned. Additionally, the parameter, 
ASSIGNMT_SOLUTION, is used as an output parameter for 
returning a solution of a lower dimensioned assignment problem 
to an instantiation of MULTI_DIM_REIiAX which is processing a 
10 higher dimensioned assignment problem. 



iy[ULTI_DIM_RELAX is initially invoked with M as the value for 
the parameter n and the PR0B_FORMULATI0N having a data 
structure (s) representing an M-dimensional assignment problem 

15 as in (2.5), in step 500 an integer m, 2<m<n is chosen such 
that the constraint sets corresponding to dimensions m+1, . . .,n 
are to be relaxed so that an m-dimensional problem formulation 
as in (2,7) results. In step 504 an initial approximation is 
provided for {u^*^, . . . ,uS} . Subsequently, in step 508 the above 

20 initial values for { u^''^ . . . , u^} are used in an iterative 
process in determining ( iP^^ . . . , iP) which maximizes 

{^mni^'"^^f ' ' ' ,u") where u^gR^^^^ and ;c=in+l, . . . ,n} for a feasible 
solution w" subject to the constraints of Problem Formulation 
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(2.7). Note that by maximizing (u'"^S . . . , u"") , some of the 
constraints that were relaxed are being forced to be satisfied 
and in so doing information is built into a solution of this 
function for solving the input assignment problem in 
"PROB_FORMULATION. " Further note that a non- smooth 

optimization technique is used here and that a preferred 
method of determining a maximum in step 508 is the bundle - 
trust method of Schramm and Zowe (H. Schramm and J. Zowe . A 
version of the bundle idea for minimizing a non- smooth 
function: Conceptual idea, convergence analysis , numerical 
results. SIAM Journal on Optimization, 2, No. 1:121-152, 
February, 1992) . This method, along with various other 
methods for determining the maximum in step 508, are discussed 
below. 

Also note that for m >2 , a solution to the optimization 
problem of step 508 is NP-hard and therefore cannot be solved 
optimally. That is, there is no known computationally 
tractable method for guaranteeing an optimal solution. Thus, 
there are two possibilities: either (a) allow m to be greater 
than 2 and use auxiliary functions similar to those disclosed 
in the first embodiment of the A:-dimensional assignment solver 
300 in section I to compute a near-optimal solution, or (b) 
make in=2 wherein algorithms such as the f orward/reverse 



94 



ft'l 

m 

Us 

Lri 



auction algorithm of D.P. Bertsekas and D.A. Castanon (D.P. 
Bertsekas and D.A. Castanon. A forward/reverse auction 
algorithm for asymmetric assignment problems. Computational 
Optimization and Applications, 1: 277-298, 1992) provides an 
5 optimal solution. 

If option (a) immediately above is chosen, then the 
auxiliary functions are used as approximation functions for 
obtaining at least a near-optimal solution to the optimization 
problem of step 508. Note that the auxiliary functions used 

10 depend on the value of m. Thus, auxiliary functions used when 
;n=3 will likely be different from those for in=4 . But in any 
case, the optimization procedure is guided by using the merit 
function, (u^/ • • •/ u"") , which can be computed exactly via a 
two-dimensional assignment problem for guiding the 

15 maximization process. 

Alternatively, if option (b) above is chosen, then two 
important advantages result, namely, the optimization problem 
of step 508 can be always solved optimally and without using 
auxiliary approximation functions. Thus, better solutions to 

20 the original M-dimensional assignment problem are likely since 
there is no guarantee that when the non-smooth optimization 
techniques are applied to the auxiliary functions the 
techniques will yield an optimal solution to step 508. 
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Furthermore, it is important to note that without auxiliary 
functions, the processing in step 508 is both conceptually 
easier to follow and more efficient. 

Subsequently, in step 512 of Fig. 8, the solution 
5 ( iP"*"^, . . . , lT") is used in determining an optimal solution to 
Problem Formulation (2.9) as generated according to (2.8) and 
Fact A. 3 . 

In step 516, the solution is used in defining the cost 
matrix c''""'''^ as in (2.12) . Subsequently, if n-jn+l=2, then the 

10 assignment problem (2.13) may be solved straightforwardly 
using known techniques such as forward/reverse auction 
algorithms- Following this, in step 528, the solution to the 
two-dimensional assignment problem is assigned to the variable 
ASSIGNMT_SOLUTION and in step 532 ASSIGNMT_SOLUTION is 

15 returned to a dimension three level recursion of 
MULTI_DIM_RELAX for solving a three-dimensional assignment 
problem. 

Alternatively, if in step 520, n--m+l>2, then in step 536 
the data structure (s) representing a problem formulation as in 
20 (2.13) is generated and assigned to the parameter, 
PROB_FORMULATION. Subsequently, in step 54 0 a recursive copy 
of MULTI_DIM_REIiAX is invoked to solve the lower dimensional 
assignment problem represented by PROB^FORMULATION. Upon the 



96 



CSUR.01USRI 

completion of step 540, the parameter, ASSIGNMT^SOLUTION, 
contains the solution to the (n-in+l) -dimensional assignment 
problem. Thus, in step 544, the (n-m+l)^^ solution is used to 
solve the n-dimensional assignment problem as discussed 
5 regarding equations (2.14) . Finally, in steps 548 and 552 the 
solution to the n-dimensional assignment problem is returned 
to the calling program so that, for example, it may be used in 
taking one or more actions such as (a) sending a warning to 
aircraft or sea facility; (b) controlling air traffic; (c) 
10 controlling anti-aircraft or anti-missile equipment; (d) 
taking evasive action; or (e) surveilling an object. 



II. 3.1. Comments on the Various Procedures 
Provided by Figure 8 

15 

There are many procedures described by Figure 8 . One 
such procedure is the first embodiment of the multidimensional 
assignment solving process of section I.l. That is, by 
specifying in=n-l in step 500, a single set of constraints is 

20 relaxed in step 508. Thus, one set of constraints is 
incorporated into the objective function via the Lagrangian 
problem formulation, resulting in an (jn=n-l) -dimensional 
problem. The relaxed problem is subsequently maximized in 
step 512 with respect to the corresponding Lagrange 

25 multipliers and then a feasible solution is reconstructed for 
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the n-dimensional problem using a two-dimensional assignment 
problem. The second procedure provided by Figure 8 is a novel 
approach which is not suggested by the first embodiment of 
section I.l. In fact, the second procedure is somewhat of a 
mirror image of the first embodiment in that n-2 sets of 
constraints are simultaneously relaxed, yielding immediately 
an (731=2) -dimensional problem in step 512. Thus, a feasible 
solution to the n-dimensional problem is then recovered using 
a recursively obtained solution to an (n-1) -dimensional 
problem via step 540. In this case, the function values and 
subgradients of (u^, . . . , u"") of step 508 can be computed 
optimally via a two-dimensional assignment problem. The 
significant advantage here is that there is no need for the 
merit or auxiliary functions, as required in the first 

embodiment of the multidimensional assignment solving process 
of section I.l above and also there is no need for the more 
general merit or auxiliary functions W^, as discussed in 
subsection II. 4. 2 below. Further, note that all function 
values and subgradients used in the nonsmooth maximization 
process are computed exactly (i.e., optimally) in this second 
procedure. Moreover, problem decomposition is now carried out 
for the n-dimensional problem; however, decomposition of the 
(n-1) -dimensional recovery problem (and all lower recovery 
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problems) is performed only after the problem is formulated. 
Between these two procedures are a host of different 
relaxation schemes based on relaxing n - m sets of constraints 
to an jn-dimensional problem (2<m<n) , but these all have the 
5 same difficulties as the procedure for the first embodiment of 
section I . 1 in that the relaxed problem is an NP-hard problem. 
To resolve this difficulty, we use an auxiliary or merit 
function as described in subsection II. 4. 2 below, (The 

notation ¥^ was used for W^.^^n section I.l.) For the case 
Cn 10 m<n-l , the recovery procedure is still based on an NP-hard 



partitioning techniques similar to those discussed in Section 
1-3 may be used to identify the assignment problem with a 
layered graph and then to find the disjoint components of this 



15 graph. In general, all relaxed problems can be decomposed 
prior to any nonsmooth computations because their structure 
stays fixed throughout the algorithm of Fig. 8. However, all 
recovery problems cannot be decomposed until they are 
formulated, as their structure changes as the solutions to the 



(n-ni+l) -dimensional assignment problem. 



Note 



that 



the 



2 0 relaxed problems change. 
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II. 4. Details and Refinements Relating to the 

Flowchart of Figure 8 

Further explanation is provided here on how various steps 

of Figure 8 are solved. Note that the refinements presented 

here can significantly increase the speed of the relaxation 

procedure of step 508. 

II. 4.1 Maximization of the Non- smooth Function 

One of the key steps in the Lagrangian relaxation 
algorithm in section II. 3 is the solution of the problem 

Maximize {^'^(u"'"^ ... ,u^) : u^eW^^^; k= m+l , . . . , n} , (2.15) 

where iio=0 for all /c=iT?+l, . . . ,n. The evaluation of 

^jnn (u""""-^, . . . / u^) requires the optimal solution of the 
corresponding minimization problem (2,7). The following 
discussion provides some properties of these functions. 

Fact A. 4. Let u'"^^, . . . , u"" be multiplier vectors associated 
with the {m+D^^ through the n^^ set of constraints (2.6), let 
0^ be as defined in (2.7), let V„(z") be the objective function 
value of the n -dimensional assignment problem in equation 
(2.6) , let z"^ be any feasible solution of (2.6) , and let be 
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an optimal solution of (2.6). Then, (u'"^^, . . . , u"") is 

piecewise affine, concave and continuous in { u""^^, . , . , u"} and 



5 (2.16) 



Furthermore, 

^n^-l^di^/U"""'/ • • wU^) < ^^(u^^^ . . .,u^) (2.17) 
for all u^gR^^^^ with uS=0 and k^m, . . . ,n. 

10 

Most of the procedures for nonsmooth optimization are 
based on generalized gradients called subgradients , given by 
the following definition . 

Definition . At u= (u"""^, . . . , u") the set d^^iu) is called a 
15 subdifferential of 0^ and defined by 

d0^{u)={q € . .xm^n^'\0^{w) -0^{u) <q^{w-u) (2.18) 

2 0 A vector q E (?0^{u) is called a subgradient . 

If z" is an optimal solution of (2.7) computed during 
evaluation of 0rnni^) f differentiating with respect to u^^ 
yields the following component of a subgradient g of 0mn(^) 
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<=E..-E E •••E<...iri (2.19) 

for i}^=lf . . . / ^A: /c=2n+I, . . . ,n. 

5 If 2"" is the unique optimal solution of (2.7), d(p^{u)={g} and 
is dif f erentiable at u. If the optimal solution is not 
unique, then there are finitely many such solutions, say 
z'^(l) , . . . , z'^(JC) . Given the corresponding subgradients, 

gr^/-../g^/ the subdif f erential (90{u) is the convex hull of 
10 {gr^,...,gr-}. 

II-4.2 Mathematical Description of 

a Merit or Auxiliary Function 

For real-time needs, one must address the fact that the 

15 non-smooth optimization problem of step 508 (Fig. 8) requires 

the solution of an NP-hard problem for m>2 . One approach to 

this problem is to use the following merit or auxiliary 

function to decide whether a function value has increased or 

decreased sufficiently in the line search or trust region 

2 0 methods: 
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^n.n("^ • • • ' • • • = (2.20) 

kn(""'"-""")' ifm = 2 

\ or {2.1) is solved optimally , 

[cD2„(u^ . . .,u'";u'""S . ..,u'') , otherwise. 



where the multipliers iP , . . . , if that appear in lower order 
5 relaxations used to construct (suboptimal ) solutions of the 
jn-dimensional relaxed problem (2.7) have been explicitly 
included. Note that ¥^ is well-defined since (2.9) can always 
be solved optimally if m=2 . For sufficiently small problems 
(2.7) or (2.9), one can more efficiently solve the NP-hard 

10 problem by branch and bound. This is the reason for the 
inclusion of the first case; otherwise, the relaxed function 
02n is used to guide the nonsmooth optimization phase. That 
the merit function provides a lower bound for the optimal 
solution follows directly from Fact A. 4 and the following 

15 fact. 



20 



Fact A. 5. Given the definition of in (2.20) above, 
¥^{iP, . . . , if;vr^\ . . . rv!^) <.0^{u''^\ . . . / U'') for all multipliers 
iP, . . .,iP,u^"^ • . .,u". (2.21) 
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The actual function value used in the optimization phase is 
however, the subgradients are computed as in (2.19), but 
with the solution ^i^.-.i^ being a suboptimal solution 
constructed from a relaxation procedure applied to the m- 
dimensional problem. Again, the use of these lower order 
relaxed problems is the reason for the inclusion of the 
multipliers {?,,..,{?". 

To explain how the merit function is used, suppose there 
is a current multiplier set (u^jd/ • • • / u^jd) it is desirable 

to update to a new multiplier set (^Se^/ * • • / ^^ew) via 
( uS^i, . . . , uLJ = ( uSld. • • • . u2id) + {Axf^', Zlu") . Then 
^mn ( i^id/ • • • / "Sid/ "Sid/ • • • / "Sid) is computed where ( 12^^^/ • • • / "Sid) is 
obtained during the relaxation process used to compute a high 
quality solution to the relaxed n? -dimensional assignment 
problem (2.7) at (if*^, . . . , u") = (u^J^, . . . , u^id) • The decision to 
accept ( "Se^r/ . • • / "ne J IS then based on 5?"™, ( "^id/ • • • / "Sid; 

uSId/ ' ■ • > "Sid) < ^mn ( "^ew/ • • • / "Sewr/ "Sew/ • • • / "new) Or SOme Other 

stopping criteria commonly used in line searches. Again, 
( u^ew/ • • • / i^eJ represents the multiplier set used in the lower 
level relaxation procedure to construct a high quality 
feasible solution to the /n-dimensional relaxed problem (2.7) 
at (u"*S...,u") = (u™;^, . . . , ujJeJ . The point is that each time 
one changes (u"*^, . . . , u") and uses the merit function 
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W^{xP , . . . f if; u^^^, , . . / u^) for comparison purposes, one must 
generally change the lower level multipliers (i?/..,,L?). 

An illustration of this merit function for jn=n-l is given 
in (A.B. Poore and N. Rijavec. Partitioning multiple data 
5 sets: multidimensional assignments and Lagrangian relaxation. 
In P.M. Pardalos and H. Wolkowicz, editors, Qqudratic 
assignment and related problems : DIMACS Series in Discrete 
Mathematics and Theoretical Computer Science, volume 16, 
pages 25-37, 1994) . 

10 

II. 4. 3 Non- smooth Optimization Methods 
By Fact A. 4 the function ^jnn(") is ^ continuous, piecewise 
affine, and concave, so that the negative of 0nm(^) is convex. 
Thus the problem of maximizing ^n^n(u) is one of nonsmooth 
15 optimization. Since there is a large literature on such 
problems, only a brief discussion of the primary classes of 
methods for solving such problems is provided. A first class 
of methods , known as subgradient methods , are reviewed and 
analyzed in (N.Z. Shor. Minimization Methods for Non- 
20 Differentiable Functions. Springer-Verlag, New York, 1985) . 
Despite their relative simplicity, subgradient methods have 
some drawbacks that make them inappropriate for the tracking 
problem. They do not guarantee a descent direction at each 
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iteration, they lack a clear stopping criterion, and exhibit 
poor convergence (less than linear) . 

A more recent and sophisticated class of methods are the 
bundle methods; excellent developments are presented in {J.-B. 
5 Hiriart-Urruty and C. Lemarechal . Convex Analysis and 
Minimization Algorithms I and JJ. Springer- Verlag, Berlin, 
1993; K.C. Kiwiel . Methods of Descent for Non-Dif f erentiable 
Optimization. In A. Dold and B. Eckmann, editors. Lecture 
Notes in Mathematics, volume 1133, Berlin, 1985. Springer- 

10 Verlag) . Bundle methods retain a set (or bundle) of 
previously computed subgradients to determine the best 
possible descent direction at the current iteration. Because 
of the nonsmoothness of 0^, the resulting direction may not 
provide a "sufficient" decrease in 0^. In this case, bundle 

15 algorithms take a "null" step, wherein the bundle is enriched 
by a subgradient close to Uj^. As a result, bundle methods are 
non-ascent methods which satisfy the relaxed descent condition 
^mn(^k+i) <^inn(^k) if ^k+i'^^k' These. methods have been shown to 
have good convergence properties. In particular, bundle 

2 0 method variants have been proven to converge in a finite 
number of steps for piecewise affine convex functionals in 
(K.C. Kiwiel. An aggregate subgradient method for non- smooth 
convex minimization. Mathematical Programming, 27:320-341, 
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1983 ; H. Schramm and J. Zowe. A version of the bundle idea for 
minimizing a non-smooth function: Conceptual idea, convergence 
analysis, numerical results. SI AM Journal on Optimization, 2, 
No. 1:121-152, February, 1992). 

All of the above non- smooth optimization methods require 
the computation of ^^(u) and a g E d0^{u) . These in turn 
require the computation of the relaxed cost coefficients 
(2.8). In both the first and second procedures discussed in 
section II. 3.1, maximization of ^2n(^) must be repeatedly 
evaluated. In the most efficient implementations presently 
known of these two procedures, it was found that at least 95% 
of the computational effort of the entire procedure is spent 
in the evaluation of the relaxed cost coefficients (2.8) as 
part of computing 02n(^) - Thus, generally a method should be 
chosen that makes as efficient use of the subgradients and 
function values as possible, even at the cost sophisticated 
manipulation of the subgradients. In evaluating three 
different bundle procedures: (a) the conjugate subgradient 
method of Wolfe used in section 1 . 1 of the first embodiment of 
the present invention; (b) the aggregate subgradient method of 
Kiwiel (K-C. Kiwiel . An aggregate subgradient method for non- 
smooth convex minimization. Mathematical Programming, 27:320- 
341, 1983) ; and (c) the bundle trust method of Schramm and 
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Zowe (H. Schramm and J. Zowe . A version of the bundle idea for 
minimizing a non-smooth function: Conceptual idea, convergence 
analysis, numerical results. SI AM Journal on Optimization, 2, 
No. 1:121-152, February, 1992), it was determined that for a 
5 fixed number of non- smooth iterations, say, ten, the bundle- 
trust method provides good quality solutions with the fewest 
number of function and subgradient evaluations of all the 
methods, and is therefore the preferred method. 



in 



10 II-4.4 The Two-Dimensional Assignment Problem 

The f orward/reverse auction algorithm of Bertsekas and 

Castafion (D.P. Bertsekas and D.A. Castafion. A f orward/reverse 

auction algorithm for asymmetric assignment problems . 

Computational Optimization and Applications, 1:277-298, 1992) 
15 is used to solve the many relaxed two-dimensional problems 

that occur in the course of execution. 



II. 4-5. Initial Multipliers and Hot Starts 
The effective use of "hot starts" is fundamental for 
20 real-time applications. A good initial set of multipliers can 
significantly reduce the number of non-smooth iterations (and 
hence the number of (P^ evaluations) required for a high 
quality recovered solution. Further, there are several ways 
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that multipliers can be reused. First, if an n -dimensional 
problem is relaxed to an m-dimensional problem, relaxation 
provides the multiplier set {u""^^, u""^^, , . . , u''} . These can be 
used as the initial multipliers for the (n-zn+l) -dimensional 
5 recovery problem for n-ni+l>2. This approach has also worked 
well to reduce the number of non- smooth iterations during 
recovery. 

Further, for track maintenance, initial feasible 
solutions are generated as follows. When a new scan of 



10 information (a new observation set) arrives from a sensor, one 
can obtain an initial primal feasible solution by matching new 
reports to existing tracks via a two-dimensional assignment 



approach. Thus, an initial primal solution exists and then we 



15 use the above relaxation procedure to make improvements to 



this TWS solution. Also for track maintenance, multipliers 
are available from a previously solved and closely related 
multidimensional assignment problems for all but the new 
observation set. 



Given a feasible solution of the multidimensional 
assignment problem, one can consider local search procedures 



problem. 



This is known as the track-while-scan (TWS) 



20 



II. 4. 6. 



Local Search Methods 
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to improve this result, as described in (C.H. Papadimitriou 
and K. Steiglitz. Comhinatorial Optimization: Algorithms and 
Complexity. Prentice-Hall, Inc., Englewood Cliffs, NJ, 1982; 
J. Pearl. Heuristics : Intelligent Search Strategies for 
5 Computer Problem Solving, Addison-Wesley , Reading, MA, 1984) . 
One method is the idea of k interchanges. Since for sparse 
problems only those active arcs that participate in the 
solution are stored, it is difficult to efficiently identify 
eligible variable exchange sets with some data structures for 

10 solving the assignment problem. However, a local search that 
may be more promising is to use the two-dimensional assignment 
solver in the following way. Given a feasible solution to the 
multidimensional assignment problem, the indices that 
correspond to active arcs in the solution are enumerated. 

15 Subsequently, one coordinate is freed to formulate a two- 
dimensional assignment problem with one index corresponding to 
the enumeration and the other to the freed coordinate, and 
then solve a two-dimensional assignment problem to improve the 
freed index position. 



II. 4. 7. Problem Decomposition 
The procedures described thus far are all based on 
relaxation. Due to the sparsity of assignment problems. 
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however, frequently a decomposition of the problem into a 
collection of disjoint components can be done wherein each of 
the components can be solved independently. Due to the setup 
costs of Lagrangian relaxation, a branch and bound procedure 
is generally more efficient for small components, say ten to 
twenty feasible 0-1 variables (i.e., ^ii-.-i^,)- Otherwise, the 
relaxation procedures described above are used. Perhaps the 
easiest way to view the decomposition method is to view the 
reports- or measurements as a layered graph. A vertex is 
associated with each observation point, and an edge is allowed 
to connect two vertices only if the two observations belong to 
at least one feasible track of observations. Given this 
graph, the decomposition problem can then be posed as that of 
identifying the connected subcomponents of a graph which can 
be accomplished by constructing a spanning forest via a depth 
first search algorithm, as described in (A.V. Aho, J.E. 
Hopcroft, and J.D. Ullman. The Design and Analysis of Computer 
Algorithms. Addison-Wesley, MA, 1974). Additionally, 
decomposition of a different type might be based on the 
identification of bi- and tri -connected components of a graph 
and enumerating on the connections. Here is a technical 
explanation. Let 

^ = {Zij,i2...iJ^iii2-"in preassigned to zero} 
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denote the set of assignable variables. Define an undirected 
graph G{N,A) where the set of nodes is 

Ar={z5jn=l, . . , ,n; . . . 

and arcs, 

5 A={ (z^^, zj^) |n^l, jn^^i Ji"^^ there exists Zi^i^,..!^^"^ 

such that jn=^n Ji=ii} . 

Note that the nodes corresponding to zero index have not 
been included in the above defined graph, since two variables 
that have only the zero index in common can be assigned 
10 independently. Connected components of the graph are then 
easily found by constructing a spanning forest via a depth 
first search. (A detailed algorithm can be found in the book 
by Aho, Hopcroft and Ullman cited above) . Furthermore, this 
procedure is used at each level in the relaxation, i.e., is 
15 applied to each assignment problem (1.4) for n=3,.--,M. 

The original relaxation problem is decomposed first. All 
relaxed assignment problems can be decomposed a priori and all 
recovery problems can be decomposed only after they are 
formulated. Hence, in the n-to- (n-l) case, we have n-2 
20 relaxed problems that can all be decomposed initially, and the 
recovery problems that are not decomposed (since they are all 
two-dimensional) . In the n-to-2 case, we have only one 
relaxed problem that can be decomposed initially. This case 
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yields r2-3 recovery problems, which can be decomposed only 
after they are formulated. 

III. NEW APPROACHES TO TRACK INITIATION AND MAINTENANCE 



III.l. Introduction 

The ever- increasing demand for sensor surveillance 
systems is to accurate track and identify objects in real- 

10 time, even for dense target scenarios and in regions of high 
track contention. The use of multiple sensors,, through more 
varied information, has the potential to greatly improve 
target state estimation and identification. However, to take 
full advantage of the available data, a multiple frame data 

15 association approach is needed. The most popular such 
approach is an enumerative technique called multiple 
hypothesis tracking (MHT) . As an enumerative technique, MHT 
can be too computationally intensive for real-time needs 
because this problem is NP-hard. A promising approach is to 

20 utilize the recently developed Lagrangian relaxation 
algorithms (K. P . Pattipati , S. Deb, and Y . Bar- Shalom. A 
multisensor-multitarget data association algorithm for 
heterogeneous sensors. IEEE Transactions on Aerospace and 
Electronic Systems, 29, No. 2:560-568, April 1993; A.B. Poore 

25 and N.Rijavec. Partitioning multiple data sets: 



5 



USING MULTIDIMENSIONAL ASSIGNMENT PROBLEMS 
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multidimensional assignments and lagrangian relaxation, in 
quadratic assignment and related problems, P.M. Pardalos and 
H.Wolkowicz, editors, DIMACS series in Discrete Mathematics 
and Theoretical Computer Science, 16:25-37, 1994; A.J. 
5 Robertson III. A class of lagrangian relaxation algorithms for 
the multidimensional problem. Ph.D. Thesis, Colorado State 
University, Ft. Collins, CO, 1995 j for the multidimensional 
assignment problem; however, there are many other potentially 
good approaches to these assignment problems such as LP 

10 relaxation combined with an interior point method, GRASP, and 
parallelization. 

These data association problems are fundamentally 
important in tracking. Thus, our first objective is . to 
present an overview of the tracking problem and the context 

15 within which these data association problems fit. Following 
a brief review of the probabilistic framework for data 
association, the second objective is to formulate two models 
for track initiation and maintenance as multidimensional 
assignment problems. This formulation is based on a window 

2 0 moving over the frames of data. The first and simpler model 
uses the same window length for track initiation and 
maintenance, while the second model uses different lengths. 
As the window moves over the frames of data , one creates a 
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sequence of these multidimensional assignment problems and 
there is an overlap region of frames common to adjacent 
windows. Thus, given the solution of one problem, one can 
develop a warm start for the solution to the next problem in 



critically important to the design of real-time algorithms. 

III. 2. Overview of the Tracking Problem 
Tracking and data fusion are complex subjects of 

10 specialization that can only be briefly summarized as they are 
related to the subject of this paper. The processing of track 
multiple targets using data from one or more sensors is 
typically partitioned into two stages or major functional 
blocks, namely, signal processing and data processing (Y.Bar- 

15 Shalom. Multitarget-Multisensor Tracking: Advanced 
Applications, Artech House, Dedham, MA, 1990; Y. Bar-Shalom and 
T.E. Fortmann, Tracking and Data Association, Academic Press, 
Boston, MA, 1988; S.S. Blackman. Multiple Target Tracking with 
Radar Applications . Artech House, Dedham, MA, 1986) . The 

20 first stage is the signal processing that takes raw sensor 
data and outputs reports. Reports are sometimes called 
observations, threshold exceedances, plots, hits, or returns, 
depending on the type of sensor. The true source of each 



5 the sequence, as shown hereinafter. 



Such information is 
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report is usually unknown and can be due to a target of 
interest, a false signal, or persistent background objects 
that can be moving in the field of view of the sensor. 

Two principal functions of data processing are tracking 
5 and discrimination or target identification. The 
discrimination or target identification function distinguishes 
between tracks that are for targets of interest and other 
tracks such as those due to background. Also, if there is 
enough information, each track is classified as to the type of 

10 target it appears to represent . Target identification and 
discrimination will not be discussed further in this paper 
except to comment here on the use of attributes (including 
identification information) in the data association. 

There are two types of attributes or features, namely, 

15 discrete valued variables and continuous valued variables. 
Available attributes of either or both types should be used in 
computing the log likelihoods or scores for data association. 
In discussing tracking in this paper, the attributes will not 
be mentioned explicitly but it should be understood that some 

2 0 reports and some tracks may include attributes information 
that should be and can be used in the track processing (both 
data association and filtering) if it is useful. 
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III. 2.1. Tracking Functions 
The tracking function accepts reports for each frame of 
data and constructs or maintains a target state estimate, the 
variance -covariance matrix of the state estimation error, and 
5 refined estimates (or probabilities) of the target attributes. 
The state estimate typically includes estimates of target 
position and velocity in three dimensions and possibly also 
acceleration and other variables of interest. 

The tracking function is typically viewed in two stages, 
10 namely, data association and filtering. Also, the process of 
constructing new tracks, called track initiation, is different 
from the process of updating existing tracks, called track 
maintenance . 

In track maintenance, the data association function 
15 decides how to relate the reports from the current frame of 
data to the previously computed tracks. In one approach, at 
most one report is assigned to each track, and in other 
approaches, weights are assigned to the pairings of reports to 
a track. After the data association, the filter function 
2 0 updates each target state estimate using the one or more (with 
weights) reports that were determined by the data association 
function, A filter commonly used for multiple target tracking 
and data fusion is the well known Kalman filter (or the 
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extended Kalman filter in the case of nonlinear problems) or 
a simplified version of it (Y. Bar-Shalom. Multitarget- 
Multisensor Tracking: Advanced Applications. Artech House, MA, 
1990; Y. Bar-Shalom and T.E. Fortmann. Tracking and Data 
Association. Academic Press, Boston, MA, 1988; S.S, Blackman. 
Multiple Target Tracking with Radar Applications . Artech 
House, Dedham, MA, 1986) . 

In track initiation, typically a sequence of reports is 
selected (one from each of a few frames of data) to construct 
a new track. In track initiation, the filtering function 
constructs a target state estimate and related information 
based on the selected sequence of reports. The new track is 
later updated by the track maintenance processing of the 
subsequent frames of reports. 

In some trackers, there is also a tentative tracking 
function for processing recently established tracks until 
there is enough confidence to include them in track 
maintenance. While for simplicity of presentation, tentative 
tracking will not be included in the processing discussed in 
this paper, the techniques discussed could readily include one 
or more tentative tracking functions. 

There have been numerous approaches developed to perform 
the data association function. Since optimal data association 
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is far too complex to implement, good but practical sub- 
optimal approaches are pursued- Data association approaches 
can be classified in a number of ways. One way to classify 
data association approaches is based on the number of data 
5 frames used in the association process (O.E. Drummond. 
Multiple sensor tracking with multiple frame, probabilistic 
data association. In Signal and Data Processing of Small 
Targets, SPIE Proceedings, volume 2561, pages 322-336, 1995) . 
In single frame data association for track maintenance, 

10 "hard decisions" are made on the assignment of reports to the 
tracks. Some single frame approaches include: individual 
nearest neighbor, PDA, JPDA, and global nearest neighbor 
(sequential most probable hypothesis tracking) , which uses a 
two-dimensional assignment algorithm (Y. Bar-Shalom. 

15 Multitarget-Multisensor Tracking: Advanced Applications , 
Artech House, MA, 1990; Y. Bar- Shalom and T.E. Fortmann. 
Tracking and Data Association. Academic Press, Boston, MA, 
1988; S.S. Blackman. Multiple Target Tracking with Radar 
Applications. Artech House, Dedham, MA, 1986) . In most single 

2 0 frame data association approaches, only one track per object 
is carried forward to be used in processing the next frame of 
reports . 
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Multiple frame data association is more complex and 
frequently involves purposely carrying forward more than one 
track per target to be used in processing the next frame of 
reports. By retaining more than one track per target, the 
5 tracking performance is improved at the cost of increased 
processing load. The best known multiple frame data 
association approach is an enumerative technique called 
multiple hypothesis tracking (MHT) which enumerates all the 
global hypothesis with various pruning rules . 

10 As shown hereinafter, the log likelihood of the 

probability of a global hypothesis can be decomposed into the 
sum of scores of each track that is contained in the 
hypothesis. As a consequence, the most probable hypothesis 
can be identified with the aid of a non- enumerative algorithm 

15 for the assignment so formulated. 



There are a number of ways that track initiation and 
track maintenance functions can interact (O.E. Drummond. 
20 Multiple target tracking with multiple frame, probabilistic 
data association. In Signal and Data Proceedings, volume 1954, 
pages 394-408, 1993) . Two ways that are especially pertinent 
to the methods of this paper will be discussed in this 
section. 



III. 2. 2. 



Interface Between Track Initiation 
And Track Maintenance 
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One approach is to treat these two functions 
sequentially, by first assigning reports to tracks and then 
using the information that remains to form new tracks. A 
better approach is to integrate both processes where 
5 initiating tracks compete for the new frame of reports on an 
equal basis, i.e., at the same time. 

In the integrated approach, the data association for 
track initiation and track maintenance processing are combined 
Q and conducted simultaneously. One assignment array is created 

10 that includes the track scores for potential tracks for both 
track initiation and track maintenance. The first dimension 
of the assignment problem includes only all the tracks either 
created or updated in the processing of the prior frames of 
reports- Each of the remaining dimensions accommodates one 
15 frame of reports. 

After the assignment algorithm finds a solution, each of 
the previously established tracks are updated with the report 
assigned to it from the second dimension of the assignment 
array. The remaining reports in the second dimension that 
2 0 were in the assignment solution are firmly established as 
track initiators. The unassigned reports in the second 
dimension of the assignment array are discarded as false 
signals. The processing of the next frame of reports repeats 
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this process using these updated tracks and the newly 
identified initiators in the first dimension of the cost array 
for processing the new frame of reports. Details of this 
approach are discussed hereinafter. 
5 The integrated approach just discussed, uses the same 

number of frames of reports for both track initiation and 
track maintenance. The goal in using the multi-dimensional 
assignment algorithm is to provide improved performance while 
minimizing the amount of processing required. Typically, 

J^; 10 track initiation will benefit from more frames of reports than 
will track maintenance. Thus, a second approach that 
integrates track initiation and track maintenance is discussed 
hereinafter, wherein the number of frames of reports is not 

^ the same for these two functions. This is a novel approach 

t\ 15 that is introduced in this paper. 

III. 2. 3. Multiple Sensor Processing 
There are many advantages and many ways of combining data 
from multiple sensors. There are also many ways of 
20 categorizing the different algorithmic architectures for 
processing data from multiple sensors. One approach outlines 
four generic algorithmic architectures (O.E. Drummond. 
Multiple sensor tracking with multiple frame, probabilistic 
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data association. In Signal and Data Processing of Small 
Targets, SPIE Proceedings, volume 2561, pages 322-336, 1995) . 
Two of these generic architectures are especially pertinent to 
this paper and are summarized briefly. 

In the Centralized Fusion algorithmic architecture, 
reports are combined from the various sensors to form global 
tracks- This algorithmic architecture is also called Central 
Level Tracking, Centralized Algorithmic Architecture, or 
simply the Type IV algorithmic architecture. In track 
maintenance, for example, if a single frame data association 
is used, then for data association a frame of reports from one 
sensor is processed with the latest global tracks; then the 
global tracks are updated by the filter function. After the 
processing of this frame of reports is completed, a frame of 
the reports from another (or the same) sensor is processed 
with these updated global tracks. This process is continued 
as new frames of data become available to the system as a 
whole . 

In Centralized Fusion, using the multi-dimensional 
assignment algorithm for the data association with multiple 
sensors is similar to the processing of data from a single 
sensor. Instead of using multiple frames of reports from a 
single sensor, the multiple frames of reports come from 
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multiple sensors. The frames of reports from all the sensors 
are ordered based on the nominal time each frame of reports is 
obtained and independent of which sensor provided the reports. 
In this way the approaches discussed in this paper can be 
extended to process reports from multiple sensors using the 
Centralized Fusion algorithmic architecture. 

The second pertinent algorithmic architecture is Track 
Fusion. This approach is also called the Hierarchical or 
Federated algorithmic architecture, sensor level tracking or 
simply the Type II algorithmic architecture. In track 
maintenance, for example, a processor for each sensor computes 
single sensor tracks; these tracks are then forwarded to the 
global tracker to compute global tracks based on data from all 
sensors . 

After the first time a sensor processor forwards tracks 
to the global tracker, then subsequent tracks for the same 
targets are cross-correlated with the existing global tracks. 
This track-to-track cross-correlation is due to the common 
history of the current sensor tracks and the tracks from the 
same sensor that were forwarded earlier to the global tracker. 
The processing must take this cross-correlation into account 
and there are a number of ways of compensating for this -cross - 
correlations. One method for dealing with this cross- 
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correlation is to decorrelate the sensor tracks that are sent 
to the global tracker. There are a variety of ways to achieve 
this decorrelation and some are summarized in recent paper 
(O.E. Drummond. Feedback in track fusion without process 
noise. In Signal and Data Processing of Small Targets, SPIE 
Proceedings, volume 2561, pages 369-383, 1995). 

Once the sensor tracks are decorrelated they can be 
processed by the global tracker in almost the same way as 
reports are processed. In the case of track fusion the 
association process is referred to as track (or track-to- 
track) association rather than data (or report-to-track) 
association. If the sensor tracks are decorrelated, the 
global tracker can process the tracks from the various sensor 
processors in much the same way that the global tracker of 
Centralized Fusion processes reports. Accordingly the methods 
described in this paper can be readily extended to processing 
data from multiple sensors using either the Centralized Fusion 
or Track Fusion or even a hybrid combination of both 
algorithmic architectures . 

III. 3. Formulation of the Data Association Problem 
The goal of this section is to briefly outline the 
probabilistic framework for the data association problems 
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presented in this work. The technical details are presented 
elsewhere (A.B. Poore . Multidimensional assignment formulation 
of data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
5 Applications, 3:27-57, 1994). The data association problems 
for multisensor and multitarget tracking considered in this 
work are generally posed (A.B. Poore. Multidimensional 
assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion, 
10 Computational Optimization and Applications, 3:27-57, 1994) as 
that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



where Z^^^ represents N+1 data sets, y is a partition of indices 
15 of the data (and thus induces a partition of the data) , r* is 
the finite collection of all such partitions, r is a discrete 
random element defined on r*, and P(r=Y|z^^^) is the posterior 
probability of a partition y being true given the data Z^^^, 
The term partition is defined below; however, this framework 
20 is currently sufficiently general to cover set packings and 
coverings . 



Maximize {p{r=y\ Z^^^) [y e r*}. 



(3.1) 
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10 

■ait.-? 



1.; 15 
: ; 



20 



Consider W+1 data sets Z{k) (k=l, . . . ,N+1) , each consisting 
of M;^ actual reports and a dummy report Zq , and let denote the 
cumulative data set defined by- 



respect ively. (The dummy report Zq serves several purposes in 
the representation of missing data, false reports, initiating 
tracks, and terminating tracks (A.B. Poore . Multidimensional 
assignment formulation of data association problems arising 
from multi target tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:21-51, 1994) .) 
In multisensor data fusion and multitarget tracking the data 
sets Z{k) may represent different classes of objects, and each 
data set can arise from different sensors. For track 
initiation the objects are reports that must be partitioned 
into tracks and false alarms. In our formulation of track 
maintenance, which uses a moving window, one data set will be 
tracks and remaining data sets will be reports which are 
assigned to existing tracks, as false reports, or to 
initiating tracks. 

We specialize the problem to the case of set partitioning 
(A.B. Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 




(3,2) 
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multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994) defined in the following way. 
Define a "track of data" as |z^y z^^^^^^j where each can 
assume zero or nonzero values. A partition of the data will 
5 refer to a collection of tracks of data wherein each report 
occurs exactly once in one of the tracks of data and such that 
all data is used up; the occurrence of dummy report is 
unrestricted. The reference partition y° is that in which all 
reports are declared to be false. 
10 Next, under appropriate independence assumptions one can 

show (A.B. Poore. Multidimensional assignment formulation of 
data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994) 

15 

^^i-^N+i ^ likelihood ratio containing probabilities for 

20 detection, maneuvers, and termination as well as probability 
density functions for report errors, track initiation and 
termination . Define 



c. .... = -InL. .... 



25 z 



^1 



J 1, if(z, .... ) are assigned to as a track, 
= \ ^1 y 1 (3.4) 



0, otherwise. 
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Then, recognizing that - In 



P(r = y|z^'M 
P(r = Y°|Z^^^) 



the problem (3.1) can be expressed as the (N+l) -dimensional 
assignment problem 



Lnimize - - - zl i i 



Min 



Subject To 



E E i,=i....,M,, 



i =0 i =0 1" 



• • • 2^ 2^ ... 2^ ~ 1 

for .,.yMj^ and k = 2,...fN, 

z 6(0,1} for all i ,..., i , 



(3.5) 



5 where Cq o is arbitrarily defined to be zero. Here, each group 
of sums in the constraints represents the fact that each non- 
dummy report occurs exactly once in a "track of data." One 
can modify this formulation to include multi -assignments of 
one, some, or all the actual reports. The assignment problem 
10 (3.5) is changed accordingly. For example, if z^^ is to be 
assigned no more than, exactly, or no less than n/^ times, then 
the "=1" in the constraint (3.5) is changed to =, ^ f " 
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respectively. In making these changes, one must pay careful 
attention to the independence assumptions, which need not be 
valid in many applications. Expressions for the likelihood 
ratios can be found in the work of (A.B. Poore. 
5 Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors. DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science^ American Mathematical Society, Providence, 
R.I., 19:169-198, 1995) and the references therein. 



111,4. Track Initiation and Maintenance 
In this section we explain two multiframe assignment 
formulations to the track initiation and maintenance problem. 
The continued use of all prior information is computationally 

15 intensive for tracking, so that a window sliding over the 
frames of reports is used as the framework for track 
maintenance and track initiation within the window. The 
objectives are to describe an improved version of a simple 
method and then to put this into a more general framework in 

2 0 which track initiation and maintenance have different length 
moving windows . 



10 
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111.4,1. The First Approach to Track 
Maintenance and Initiation 



The first method as explained in this section is an 
improved version of our first track maintenance scheme (A.B. 
5 Poore . Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994) and uses the same window length 
for track initiation and maintenance after the initialization 

10 step. The process is to start with a window of length i^+1 
anchored at frame one. In the first step there is only one 
track initiation in that we assume no prior existing tracks. 
In the second and all subsequent frames, there is a window of 
length N anchored at frame k plus a collection of tracks up to 

15 frame k. This window is denoted by {k, k-hl , . . . , k-^N) . The 
following explanation of the steps is much like mathematical 
induction in that we explain the first step and then step k to 
step k+1. 



be an enumeration of all those zero-one variables in the 
solution of the assignment problem (3,5) (i.e., z. . . =1 ) 



excluding all the false reports in the solution (i.e., all 



Track Maintenance and Initiation: Step 1. Let 



20 




(3.6) 
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those zero-one variables with exactly one nonzero index) and 
zero-one variables in the solution for which (ii,i2 )-(0,0}. 
(The latter can correspond to tracks that initiate on frames 
three and higher.) These denote our initial tracks. 
5 Consider only the first two index sets in this 

enumeration (3.6) and add the zero index l2=0 with the 
corresponding values of i^ and ±2 being zero. Thus, the 

enumeration is now {-^1 (-^2^ ' ^2 ^-^2) }i'=o ' "^^^ notation 

r_(I,) = (z. .,z. J will be used for this pairing. Suppose 
10 now that the next data set, i.e., the set, is added to 

the problem. 

To explain the costs for the new problem, one starts with 
the hypothesis that a partition y^r* being true is now 
conditioned on the truth of the pairings on the first two 
15 frames being correct. Corresponding to the sequence 
jr^ , z. , z. f, the likelihood function is then given by 

L = n L, . . , Inhere 

Y , . ■^2^3-" ^N*2 



{r,(i,),r,^,...,z,^^^}eY (3.7) 



L„ ,„ =1. 



Next, define the cost and the corresponding zero-one variable 



20 
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C, . . = -InL, . . , 

(3.8) 



Jl, if\zl,...,z^;;fj is assi 
^2V^».2 otherwise, 



gned to Track 72 (Ij) ^ 



respectively. Then the track maintenance problem Maximize 
{Lj,|Y£^r*} can be formulated as the following multidimensional 
assignment 

10 Minimize 22 z2 • ■ • z2 i i ± i 

Subject to: ^ tkt^ 



(3.9) 



i3=0 i„,2 = 0 



• • • zJ ^7i...i =1' 23 = 1,. ..,M, 



1^=0 ±3=0 ip.i = 0 ip,i = 0 1^,2 = 0 



for ip=l,...,Mp and p=4,..., W+2-1, 



V3..-i,,,^{0'l} for all I2, 23, -,1^,2 



Track Maintenance and Initiation: Step k. At the 

beginning of the k^^ step, we solve the following (N+1) - 
15 dimensional assignment problem. 
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Minimize ^ • • • ^v.r-i./i.i.. ' 



1„=0 i,.i = 0 i,,„=0 



Subject to: 



2^ 2^ ' ' • 2^ 2^ ' ' • 2^ ^ 1 i i ^ 



1,=0 i,,i = o ip.i = o ip,i = 0 i,,^=0 



for ip=l,.„,M^ and p=k+2, ...,N-^k-l, 



Li, M. 



^ ^^k^-2 

L =0 i =0 i =0 * 

^k ^ -^k+l ^ •^k+l+N-2 ^ 



where for the first step li and are replaced by and M^, 
respectively. The first index 1^ in the subscripts corresponds 
to the sequence of tracks {Tj^ ( ij;^) |^^_^, where 
r^(l^) =|z/^ (1^) , zj^^(ljj^) I is a track of data from the solution 
10 of the problem prior to the formulation of (3.10). If k=l, 
then it is just the first data set in frame one. 

Basic Assxjmption. Suppose problem (3.10) has been solved 
and let the solution, i.e., those zero-one variables equal to 
one, be enumerated by 

15 ((i. ' (^..i) ' )}i';;=i (3.11) 

with the following exceptions. 
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a. All zero-one variables for which (1^, i^^^^) = (0, 0) are 
excluded. Thus, tracks that initiate on frames after the 
(k+1)^^ are not included in the list. 

b. All zero-one variables whose subscripts have the form 1^=0 
and exactly one nonzero index in the remaining indices 
{i^^^, i^^^) are excluded. These correspond to false 
reports on frames p = k-\-l , , k+N , 

c. All variables z, , . for which i., J = (0, 0,..., 0) 
and ij^(I^)=0 are excluded from (1.11). In other words, 
the reports on the last N+1 frames in the 
[Tj^{1j^) , z^^^, . . . , z^^^^ are all dummy. Any solution with 
this feature corresponds to a terminated track. 

Given the enumeration (3.11), one now fixes the assignments 
on the first two index sets in the list (3.11). The zero 
index Ij^+i'^ is added to the enumeration to specify 
(i^(0) , i^^^ (0) )= (0, 0) that is used to represent false reports and 
tracks that initiate on frame k+2 or later, so that the 

enumeration (3.11) is now {(-^/c ^ -^/c+i^ ' -^/c+i ^ ^;c+iO}i|^'|=o * 

Then, for 1^ - =1, L. . , the l/^ such track is denoted by 

X+i >C+1 fc+1 

r..i(Vi)= <;(^.n)} an<i the (N+l)-tuple 

{r..i(-^.n)'<''-'<LT} ^^11 ^®^°te a track T,^,{1,^,) plus a set 
of reports l^iVa ' '"' ^■CuTI' actual or dummy, that are feasible 
with the track T^^^il^^O ■ The (W+1) -tuple [Tj^^JO) , z^*l, z.*]^"] 
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will denote a track that initiates in the sliding window, 
i.e., on subsequent frames. A false report in the sliding 
window is one with all but one non-zero index ip for some 
p=k+2,.„,k+I+Ar in the (N+1) -tuple [t^^^{Q) , zf^^^, z^^^^^^^^^ 

The corresponding hypothesis about a partition y E V 
being true is now conditioned on the truth of the L^+i tracks 
existing at the beginning of the N-frame window. (Thus the 
assignments prior to this sliding window are fixed.) The 
likelihood function is given by 



L- n L, . , ivhere 



Next, define the cost and the zero-one variable by 



h i ...i = -InL = -InL n )z z 



. < 1' {<''•••' <LT} assigned to T^^.H,^,) , 



(3.13) 



0, otherwise, 



respectively, so that the track extension problem, which was 
originally formulated as Maximize [l Iy^T*}, can be expressed 
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as exactly the same multi-dimensional assignment in (3.4) with 
k replaced by Jc+l . Thus, we do not repeat it here- 
in. 4. 2. The Second Approach to Track 
5 Maintenance and Initiation 

In the first approach the window lengths were the same 
for both track maintenance and initiation. This can be 
inefficient in that one can usually use a shorter window for 
track maintenance than for track initiation. This section 

10 addresses such a formulation. 

The General Step To formulate a problem for track 

initiation and maintenance we consider a moving window 
centered as frame k of length J+J+1 denoted by 
[Jc-I, Jc, A:+J"] . In this representation, the window length 

15 for track maintenance is J and that for initiation, J+J+l. The 
objective will be to explain the situation at this center and 
then the move to the same length window at center k+l denoted 
by [^+1 -T, ... , Jc+l , ... , Ar+l+J"] , i.e., by moving the window to the 
right one frame at time. The explanation from the first step 

20 follows hereinafter. 

The notation for a track of data is 
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where the index Ij^ is used for an enumeration of those reports 
paired together. We also use the notation "^p^k^-^k^ to denote 
the sequence of reports belonging to track ^jt^-^jt^ but 
restricted to frames prior to and including p. Thus, 

Tp.k^h) (3.15) 

Given this notation for the tracks and partition of the 
data in the frames {k- 1, ...^ k-^J} , ^rp,jt(ijt) will denote the 
accumulated likelihood ratio up to and including frame p(p<k) 
for a track that is declared as existing on frame as a 
10 solution of the assignment problem. In this notation, the 
likelihood for '^p,k^-^k'^ ^i^d that of the association of 

I z/^^^^, z^^"^^ I with track '^k^-^k'^ is given by 



ii ) "^r (i (I (i ) any p ^ k-1, 

(2 )i ...i (i)^i (1 ) ..z (I )i ...i for any k-1, 

^k^-^k'-^kfl •^k*N ^p,k^'^k' "'■p+l^-^Jfe' ^Jt ^ -^Jt' "^Jt+W 



15 respectively. 
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The cost for the assignment of | z/^"^^, --v ^^^^"^^ ^ to track 
T^[l^) and the corresponding zero-one variable are given by 

^k^k^i ^k*j ^k^-^k* ^k^i"'^k^j and 

" '-^^^^k^h) ^^k) "^k<h) ^k.i-^kJ 

(3.17) 

1, if \ z^^'^ , zf""^ \ is assigned to track T. (i. ) , 
0, otherwise, 



respectively. Likewise, for costs associated with the false 
reports on the frames k-I to k and as associated with 
10 I z^?^^"^^ , ... , z/^"^^ ^ and the corresponding zero-one variables are 
given by: 



15 



- -ln[L, 



1, if { z z z } 

■^k-I -^k-^k+X -^k+J 

are assigned as a track, 
0, otherwise. 



] 



(3.18) 



For the window of frames in the range {k-I,.,.,k}, the easiest 
explanation of the partition of the reports is based on the 
definitions 
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p, k 



i I zf belongs to one of the tracks 



listed on frame k, i.e., 
of the T,{1,) for some 1, 



to one 
"If • • ' f j^-' 



(3.19) 



so that all of those reports not used in the tracks '^k^-^k^ on 
frame p are put into the set of available reports ^p^j^ for the 
formation of tracks over the entire window. Thus, we 
formulate the assignment problem as 



Minimize 22 52 z2 zl i i ± ^1 i i i 

^ . '^T^i ' ^k-r"-^k'^k*\--^k*J ■^k-r"-^k-^k*l--^k 

p=k-j ipefp^jt r=k^l 1^=0 

^h^k*i-^k*/h^k^i-^k*j 



1^=1 r=Jt+l i,= 0 



Subject To i E S^f ii i =1 (3.20) 

for i^GF^^j\{0] and q^k-I,.,.,k, 

for i^=l^...,M^ and q= k+1, k-^J, 
z 6{0,1} for all i i , i i , 

■^k-r"^k-^k*V"'^k^J ^ K K^i K^u 

W ^ J, 



for q=k+l, k+J, 
z . . e {0,1} for all l^^l, i^^^,..., i^ . 
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Basic Assumptions. Suppose problem (1.20) has been solved 
let the solution, i.e. those zero-one variables equal to 
, be enumerated by 



L the following exceptions. 
All zero-one variables in the second list for which 

are excluded. Thus, tracks that 
initiate on frames after the (Jc+l) are not included in 
the list. 

All false reports are excluded, i.e., all zero-one 
variables in the second list whose subscripts have 
exactly one nonzero index. 

All variables ^ ^ for which (i^+if '"^ ^jt+J 0' 0) 
and Z{p) r\T^{l^) ^0 for p^k-K,...,k where K >0 is user 
specified. Thus the track "^k^^k^ is terminated if it is 
not observed over iC+J+l frames, 
in the enumeration (3.20), one now fixes the assignments on 
all index sets up to and including the index sets. 



5 




(3.21) 
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T (1 ) =' 



f k-I k k+1 I 



(3.22) 



Thus one can now formulate the assignment problem for the 
next problem exactly as in (3.20) but with k replaced by k+1. 
Thus, we do not repeat it here. 

fi 

The Initial Step. Here is one explanation for the 

ri 

5 initial step. First, assume that N=I-^J. In this case, we 

f\ 

y Start the track initiation with a solution of (3.5) . Let 

U be an enumeration of the solution set of (3,5), i.e., those 

J zero-one variables , including Zoo...o=l 

10 corresponding to 2^=0 , but excluding all those zero -one 
variables that are assigned to one and correspond to false 
reports (i.e., there is exactly one nonzero index in the 
subscript of ^i^i^-i^^^) , all those zero-one variables that are 
assigned to one and correspond to tracks that initiate on 
15 frames higher than J+2. Then we fix the data association 
decisions corresponding to the reports in our list of tracks 
prior to and including frame Jc+I=J+2, This defines the k for 
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problem (3.5) and one can then continue the development by 
adding a frame to the window as in the general case. 

If J+J" > N, then one possibility is to start the process 
with N+1 frames, and assuming J<N, proceed as before replacing 

J by N-J for the moment, and continue to add frames without 
lopping off the first frame in the window until reaches a 
window of length I+J+l. Then we proceed as in the previous 
paragraph . 

If J+J < N, then one can solve the track initiation 
problem (3.5), formulate the problem with the center of the 
window at ic+I = N-^l-J, enumerate the solutions as above, and 
lop off the first N-J-I frames. Then, we proceed just as in 
the case I+J=N. 

III. 5. CONCLUDING COMMENTS 

A primary objective in this work has been to demonstrate 
how multidimensional assignment problems arise in the tracking 
environment. The problem of track initiation and maintenance 
has been formulated within the framework of a moving window 
over the frames of data. The solution of these NP-hard, 
noisy, large scale, and sparse problems to the noise level in 
the problem is fundamental to superior track estimation and 
identification. Thus, one must utilize the special structure 
in the problems as well take advantage of special information 
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that is available. Since these moving windows are 

overlapping, there are some algorithm efficiencies that can be 
identified and that take advantages of the overlap in the 
windows from one frame of reports to the next. Here is an 
example of the use of a primal solution of one problem to warm 
start the solution of the next problem in the sequence. 

Suppose we have solved problem (3.10) and have enumerated 
all those zero-one variables in the solution of (3.10) as in 
(3.11). Add the zero index Ij^^j^^Q ^ so that the enumeration is 

With this enumeration one can define the cost by 

^jt+i-^Jt+i+w -^Jt^-^ic+i^ *^Jt+i^ -^w+i ^-^Jt+i^ ^Jt+l+W 

and the two-dimensional assignment problem 

Minimize 22 z2 i i ^V^{z^} 

Subject to 22 ^k^i^^'-'^k^i' 

V* 2 =1 • =1 M (3.25) 

■'■k'H'^ 

for all 

Let w be an optimal or feasible solution to this two- 
dimensional assignment problem and define 
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^i.nW. = ^ ^^"^^ = (3.26) 

or if (i,,,,i,^,^„) = {0,0) 

0, otherwise. 

This need not satisfy the constraints in that there are 
5 usually many objects left unassigned. Thus, one can complete 
the assignment by using the zero-one variables in (3.10) with 
k replaced by k+l with exactly one nonzero index corresponding 
to any unassigned object or data report. 

For the dual solutions, the multipliers arising from the 

10 solution of the two-dimensional assignment problem (3.25) 
corresponding to the second variable, i.e., 1 ^ij^^i.^ I 
These are good initial values to use in a relaxation scheme 
(A.B. Poore and N.Rijavec. Partitioning multiple data sets:, 
multidimensional assignments and Lagrangian relaxation, in 

15 quadratic assignment and related problems; P.M. Pardalos and 
H. Wolkowicz, editors, DIMACS series in Discrete Mathematics 
and Theoretical Computer Science, 16:25-37, 1994; A.J. 
Robertson III. A class of lagrangian relaxation algorithms for 
the multidimensional assignment problem. Ph.D. Thesis, 

2 0 Colorado State University, Fr. Collins, CO, 1995) . Finally, 
note that one can also develop a warm start for problem (3.20) 
in a similar fashion. 
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IV. REVIEW OF NEW RELAXATION SCHEMES 



I V . 1 . Introduction 



Multidimensional assignment problems govern the central 
problem of data association in multisensor and multitarget 
5 tracking, i.e., the problem of partitioning observations from 
multiple scans of the surveillance region and from either 
single or multiple sensors into tracks and false alarms. This 
fundamental problem can be stated as (A.B. Poore. 
Multidimensional assignments and multitarget tracking, 



10 partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 



editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I., 19:169-198, 1995; A.B. Poore . Multidimensional 
assignment formulation of data association problems arising 



15 from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994) 



20 



25 
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Minimize 



t...t 




Subject to: 




(4.1) 



2^ • • • 2-^ • • • 2^ 

i, = 0 ip., = 0 ip,, = 0 i^=0 




^ • • • 2^ 




where Cq-.-o is arbitrarily defined to be zero and is included 
for notational convenience. One can modify this formulation 
to include multi -assignments of one, some, or all of the 



missing data, false alarms, initiating and terminating tracks 
as described in (A.B. Poore . Multidimensional assignments and 
multitarget tracking, partitioning data sets, I.J. Cox, P. 
Hansen, and B. Julesz, editors, DIMACS Series in Discrete 
Mathematics and Theoretical Computer Science ^ American 
Mathematical Society, Providence, R.I., 19:169-198, 1995; A.B. 
Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994). In these problems, we assume 



actual reports . 



The zero index is used in representing 
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that all zero-one variables z . . with precisely one nonzero 
index are free to be assigned and that the corresponding cost 
coefficients are well-defined. (This is a valid assumption in 
the tracking environment (A.B. Poore . Multidimensional 
5 assignments and multitarget tracking, partitioning data sets, 
I.J. Cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R.I,, 19:169-198, 
■'^^ 1995; A.B. Poore. Multidimensional assignment formulation of 

•.I ' 

10 data association problems arising from multitarget tracking 

''ll and multisensor data fusion. Computational Optimization and 

''l^ Applications, 3:27-57, 1994).) Although not required, these 

L:; cost coefficients with exactly one nonzero index can be 

ifl 

translated to zero by cost shifting (A.B. Poore and N. 

: an 
: ' 3 

15 Rijavec. A lagrangian relaxation algorithm for 
multidimensional assignment problems arising form multitarget 
tracking. SIAM Journal of Optimization, 3, No. 3:544-563, 
1993) without changing the optimal assignment. Finally, this 
formulation is of sufficient generality to include the 

20 symmetric problem and the asymmetric inequality problem (A.J. 
Robertson III. A class of lagrangian algorithms for the 
multidimensional assignment problem. Ph.D. Thesis, Colorado 
State University, Fr. Collins, CO, 1995) . 
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The data association problems in tracking that are 
formulated as Equation (4.1) have several characteristics. 
They are normally sparse, the cost coefficients are noisy and 
the problem is NP-hard (M.R. Galey and D.S. Johnson. Computers 
and Intractability. W.H. Freeman and Company, San Francisco, 
CA, 1979), but it must be solved in real-time. The only known 
methods for solving this NP-hard problem optimally are 
enumerative in nature, with branch-and-bound being the most 
efficient; however, such methods are much too slow for real- 
time applications. Thus one must resort to suboptimal 
approaches. Ideally, any such algorithm should solve the 
problem to within the noise level, assuming, of course, that 
one can measure this noise level in the physical problem and 
that the mathematical method provides a way to decide if the 
criterion has been met. 

There are many algorithms that can be used to construct 
suboptimal solutions to NP-hard combinatorial optimization 
problems. These include greedy (and its many variants), 
relaxation, simulated annealing, tabu search, genetic 
algorithms, and neural network algorithms (C.H. Papadimitriou 
and K.Steiglitz. Combinatorial Optimization: Algorithm and 
Complexity. Prentice-Hall, Inc., Englewood Cliffs, NJ, 1982; 
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J. Pearl. Heuristics: Intelligent Search Strategies for 
Computer Problem Solving. Addison-Wesley, Reading, MA, 1984; 
C.R. Reeves ed. Modern Heuristic Techniques for Combinatorial 
Problems. Halstead Press, Wiley, NY, 1993) , For the three- 
dimensional assignment problem, Pierskalla (W. Pierskalla. The 
tri-substitution method for three-dimensional assignment 
problem. Journal du CORS, 5:71-81, 1967) developed the tri- 
substitution method, which is a variant of the simplex method. 
Frieze and Yadegar (A.M. Frieze and J. Yadegar. An algorithm 
for solving 3 -dimensional assignment problems with application 
to scheduling a teaching practice. Journal of the Operational 
Research Society, 39:989-955, 1981) introduced a method based 
on Lagrangian relaxation in which a feasible solution is 
recovered using information provided by the relaxed solution. 
Our choice of approaches is strongly influenced by the need to 
balance real-time performance and solution quality . 
Lagrangian relaxation based methods have been used 
successfully in prior tracking applications (K.R. Pattipati, 
S. Deb, and Y. Bar-Shalom. A s-dimensional assignment 
algorithm for track initiation. In Proceedings of the IEEE 
Systems Conference, Kobe, Japan, pages 127-130, 1992; Y. Bar- 
Shalom, S. Deb, K.R. Pattipati, and H. Tsanakis. A new 
algorithm for the generalized multidimensional assignment 
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problem. In Proceedings of the IEEE International Conference 
on Systems, Math, and Cybernetics, Chicago, pages 132-136, 
1992; A.B. Poore and N. Rijavec. A numerical study of some 
data association problems arising in multitarget tracking. 
5 Large Scale Optimization: State of the Art, W.W. Hager, D,W, 
Hearn and P.M. Pardalos, editors. Kluwer Academic Publishers 
B.V., Boston, pages 339-361, 1994; A.B. Poore and N. Rijavec. 
Partitioning multiple data sets: multidimensional assignments 
and lagrangian relaxation. In P.M. Pardalos and H. Wolkowicz, 

10 editors , Quadratic assignment and related problems : DIMACS 
Series in Discrete Mathematics and Theoretical Computer 
Science, volume 16, pages 25-37, 1994; A.B. Poore and N. 
Rijavec. A lagrangian relaxation algorithm for 
multidimensional assignment problems from multitarget 

15 tracking. SIAM Journal of Optimization, 3 , No. 3:544-563, 
1993) . An advantage of these methods is that they provide 
both an upper and lower bound on the optimal solution, which 
can then be used to measure the solution quality. These works 
extend the method of Frieze and Yadegar (A.M. Frieze and 

2 0 J. Yadegar. An algorithm for solving 3 -dimensional assignment 
problems with application to scheduling a teaching practice. 
Journal of the Operational Research Society, 32:989-995, 1981) 
to the multidimensional case. 
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IV. 2. Probabilistic Framework for Data Association. 

The goal of this section is to explain the formulation of 
the data association problems that governs large classes of 
data association problems in centralized or hybrid 
centralized- sensor level mult i sensor /multitarget tracking . 
The presentation is brief; technical details are presented for 
both track initiation and maintenance in the work of (A.B. 
Poore. Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I., 19:169-198, 1995; A.B. Poore. Multidimensional 
assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994). 
The formulation is of sufficient generality to cover the MHT 
work of Reid, Blackman and Stein (S.S. Blackman. Multiple 
Target Tracking with Radar Applications . Artech House, 
Norwood, MA, 1986) and modifications by Kurien (T.Kurien. 
Issues in the designing of practical multitarget tracking 
algorithms. In Multitarget -Multi sensor Tracking: Advanced 
Applications by Y. Bar-Shalom. Artech House, MA, 1990) to 
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include maneuvering targets. As suggested by Blackman (S.S. 
Blackman. Multiple Target Tracking with Radar Applications. 
Artech House, Norwood, MA, 1986) , this formulation can also be 
modified to include target features (e.g., size and type) into 
the scoring function. The recent work (A.B. Poore and O.E. 
Drummond. Track initiation and maintenance using 
multidimensional assignment problems. In D.W. Hearn, W.W. 
Hager, and P.M. Pardalos, editors. Network Optimization, 
volume 450, pages 407-422, Boston, 1996. Kluwer Academic 
Publishers B.V.) significantly extends the work of this 
section to new approaches for multiple sensor centralized 
tracking. Future work will involve extensions to track-to- 
track correlation. 

The data association problems for multisensor and 
multitarget tracking considered in this work are generally 
posed as that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



where represents N data sets, y is a partition of indices 
of the data (and thus induces a partition of the data) , r* is 
the finite collection of all such partitions, r is a discrete 



Maximize 




(4.2) 
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random element defined on r*, and P(r = yI^^) is the posterior 
probability of a partition y being true given the data Z^. The 
term partition is defined below; however, this framework is 
currently sufficiently general to cover set packings and 
coverings (A.B. Poore . Multidimensional assignment formulation 
of data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994; A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R.I,, 19:169-198, 
1995) . 

Consider N data sets Z{k) (Jc=l,...,N) each of Mj, reports 
l^ij^, J ^/ and let denote the cumulative data set defined by 

^^^^={^i,f\ and Z^={Z(1),...,Z(W)}, (4.3) 

respectively. In multisensor data fusion and multitarget 
tracking the data sets Z{k) may represent different classes of 
objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must be 
partitioned into tracks and false alarms. In our formulation 
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of track maintenance (A.B. Poore. Multidimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994; 
A.B- Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J. cox, P. Hansen, and B. 
Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, American Mathematical Society, 
Providence, R.I., 19:169-198, 1995), which uses a moving 
window over time, one data set will be tracks and remaining 
data sets will be measurements which are assigned to existing 
tracks, as false measurements, or to initiating tracks. In 
sensor level tracking, the objects to be fused are tracks 
(S.S, Blackman. Multiple Target Tracking with Radar 
Applications . Artech House, Norwood, MA, 1986) . In 
centralized fusion (S.S. Blackman. Multiple Target Tracking 
with Radar Applications , Artech House, Norwood, MA, 1986) , the 
objects may all be measurements that represent targets or 
false reports, and the problem is to determine which 
measurements emanate from a common source. 

We specialize the problem to the case of set partitioning 
(A.B. Poore. Multidimensional assignment formulation of data 
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association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
Applications f 3:27-57, 1994; A.B. Poore . Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
5 I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R,I., 19:169-198, 
1995) defined in the following way. First, for notational 
convenience in representing tracks, we add a zero index to 



10 each of the index sets in Equation (4.3) , a dummy report Zq to 
each of the data sets Z{k) in Equation (4.3), and define a 



value of 0 . A partition of the data will refer to a collection 
of tracks of data wherein each report occurs exactly once in 



15 one of the tracks of data and such that all data is used up; 
the occurrence of dummy report is unrestricted. The dummy 
report serves several purposes in the representation of 

missing data, false reports, initiating tracks, and 
terminating tracks (A.B. Poore. Multidimensional assignment 



20 formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994; 
A.B. Poore. Multidimensional assignments and multitarget 



"track of data" as 




where can now assume the 
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tracking, partitioning data sets, I.J. cox, P. Hansen, and B. 
Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, American Mathematical Society, 
Providence, R.I., 19:169-198, 1995). The reference partition 
5 is that in which all reports are declared to be false. 

Next under appropriate independence assumptions one can 
show (A.B. Poore. Multidimensional assignment formulation of 
data association problems arising from multitarget tracking 

10 and multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994; A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer . Science, 

15 American Mathematical Society, Providence, R.I,, 19:169-198, 
1995) 



where y° is a reference partition, is a likelihood ratio 

20 containing probabilities for detection, maneuvers, and 
termination as well as probability density functions for 
measurement errors, track initiation and termination. Then if 



p(r = Y|z ^) = ^ n 
p(r = Y'|z^) ^ ^ ^ (ir-iJ^Y 




(4.4) 



c, 



= -InL 



N 
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-In 



p(r = Y|z ^) 
p(r = Y°|z^) 



(i,,...,i^)6Y ^ ^ 



5 Using Equation (4.4) and the zero-one variable 

^i^-i^ = 1 if (i^,-,i^) 6 Y and 0 otherwise, one can then write 

the problem (4.4) as the following i\7'-dimensional assignment 
problem: 

Minimize X! c. . z. . 

Subject To 52 i =1 {L=l,...rM), (4.6) 
E i =1 (i„=l,...,M^ a/3d p = 2,...,N-l) , 

. E (i„=i,...,M„), 

' i i IN 



10 where Cq.„q is arbitrarily defined to be zero. Here, each group 
of sums in the constraints represents the fact that each-non- 
dummy report occurs exactly once in a "track of data." One 
can modify this formulation to include multi -assignments of 
one, some, or all the actual reports. The assignment problem 

15 Equation (4.6) is changed accordingly. For example, if is 
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to be assigned no more than, exactly, or no less than 
times, then the " = 1" in the constraint (4.6) is changed to 
<, =, ^k^i^," respectively. Modifications for group tracking 

and mult i -resolution features of the surveillance region will 
5 be addressed in future work. In making these changes, one 
must pay careful attention to the independence assumptions 
that need not be valid in many applications. 

Expressions for the likelihood ratios can be found 

in the work of Poore (A.B, Poore . Multidimensional assignment 

10 formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994; 
A.B. Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J. cox, P. Hansen, and B. 

15 Julesz, editors, DIMACS Series in Discrete Mathematics and 
Theoretical Computer Science, American Mathematical Society, 
Providence, R.I., 19:169-198, 1995). These expressions 
include the developments of Reid, Kurien (T. Kurien. Issues in 
the designing of practical multitarget tracking algorithms. In 

2 0 Multitarget -Multi sensor Tracking: Advanced Applications by Y. 
Bar-Shalom, Artech House, MA, 1990) , and Stein and Blackman 
(S.S. Blackman. Multiple Target Tracking with Radar 
Applications . Artech House, Norwood, MA, 1986) . What's more, 
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they are easily modified to include target features and to 
account for different sensor types. In track initiation, the 
i\r data sets all represent reports from N sensors, possibly all 
the same. For track maintenance, we use a sliding window of 
5 N data sets and one data set containing established tracks. 
(A.B. Poore . Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994; A.B. Poore. Multidimensional 

10 assignments and multitarget tracking, partitioning data sets, 
I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R.I., 19:169-198, 
1995) . The formulation is the same as above except that the 

15 dimension of the assignment problem is now N+l . 

IV. 3. Overview of the Lagrangian Relaxation Algorithm, 

Having discussed the iV-dimensional assignment problem 
(3.1), we now turn to a description of the Lagrangian 
20 relaxation algorithm. The algorithm will proceed iteratively 
for a loop k=l, , . ., N-2 . At the completion, there remains 
one two-dimensional assignment problem that provides the last 
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Step which yields an optimal (sometimes) or near-optimal 
solution to the original N-dimensional assignment problem. In 
step k of this loop (summarized in Section IV. 3) one starts 
with the following {N-k+l) -dimensional assignment problem with 
one change in notation. If k=l, then the index notation I-^ and 
L2 are to be replaced by and M^, respectively. 

... , . w-Jt+l w-Jt+i 

Minimize 2^ ^ i i i 

Subject To J3 ^iT^ i =1, 1=1, . . . ,L., (4.7) 

E N-kU =1 i =1 M 

for i =1,...,M and p^k-^2, . . . ,N-1, 

. . ■^V*.i---i« ^' • • • '^V 

z 6 {0,1} for all 1 , i , . . . , i , 

To ensure that a feasible solution of Equation (4.7) always 
exists for a sparse problem, all variables with exactly one 
nonzero index (i.e., variables of the form ^I'o^^^o for 
1}^=1, . . . fL^ and z Qli^Q.„o for ip=I/ . . . and p=Jc+I, . . . ,N) are 

assumed free to be assigned with the corresponding cost 
coefficients being well-defined. This assumption is valid in 



CSUR.01USRI 

the tracking environment (A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
American Mathematical Society, Providence, R.I,, 19:169-198, 
1995; A.B. Poore, Multidimensional assignment formulation of 
data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994). 

Section IV. 3.1 presents some of the properties associated 
with the Lagrangian relaxation of (4.7) based on relaxing the 
last (N-k) -sets of constraints to a two-dimensional one. 
Section IV. 3. 2 describes a new approach to the problem of 
recovering a high quality feasible solution of the original 
(i\7'-Jc+l ) -dimensional problem given a feasible solution (optimal 
or suboptimal) of the relaxed two-dimensional problem is 
described hereinafter. A summary of the relaxation algorithm 
is given in Section IV. 3. 3, and in Section IV. 3. 4, we 
establish the maximization of the Lagrangian dual (an 
important aspect of the relaxation procedure) to be an 
unconstrained nonsmooth optimization problem and then present 
a method for computing the subgradients . 
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IV. 3.1 The Laqranqian Relaxed Assignment Problem 
The (N-Jc+l) -dimensional problem (4.7) has [N-k-^-l) sets 
of constraints. A (Mp+1) -dimensional multiplier vector (i.e., 
vP E R^p'^'^) associated with the p^^ constraint set will be 
5 denoted by u^= (u^, uf, . , . , ) with = 0 being fixed for each 
p = k 2,...,2\r and included for notational convenience only. 
Now, the (N-k+l) -dimensional assignment problem (4.7) is 
relaxed to a two-dimensional assignment problem by 
incorporating {N-k-1) sets of constraints into the objective 
10 function via the Lagrangian. Although any (N-k-l) constraint 
sets can be relaxed, we choose the last {N-k-1) sets of 
constraints for convenience. The relaxed problem is 
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<^N-jc.i("''''' ^ Minimize ( z"-*^^• u*^^^ . . . u ") 

... . . w-Jt+1 w-;c+i 

Minimize 7, c, , 2? 



p=;c+2 ip=o 



t "id E 



Minimize 

- £ t 



^1 i ... i 



p=;c-2 ip=o 



Subject To 53 ^Ti^'^i -i = 1/ -^^=1' - • • ^j^/ 



W-Jt + 1 



(4.8 



One of the major steps in the algorithm is the maximization of 
( u^""^, . . . , u ^) with respect to the multipliers 
(u^^^, . . . , li . It turns out that ^^-jt+i is a concave, 
continuous, and piecewise affine function of the 
5 multipliers (u^"^, . . . , u ^) , so that the maximization of ^^^^^^ is 
a problem of nonsmooth optimization. Since many of these 
algorithms require a function value and a subgradient of ^^.j^^^ 
at any required multiplier value (u^^^, . . u ^) , we address 
this problem in the next subsection. We note, however, that 
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there are other ways to maximize ^nd the next subsection 

just addresses one such method. 



up 



m 
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IV. 3. 2 Properties of the Lagrangian Relaxed Assignment Problem 
5 For a function evaluation of ^^^.j^+i / show that an 

optimal (or suboptimal) solution of this relaxed problem (4.8) 
can be constructed from that of a two-dimensional assignment 
problem. Then, the nonsmooth characteristics of ^^.j^^i 
addressed, followed by a method for computing the function 
10 value and a subgradient . 

Evaluation of ^^^j^^i' Define for each (I;,, i;,^^) an index 
(jk^2f-fjN) = (jk^2(lkf ik^i) /•"/ jw(i;c/i;c+i) ) and a new cost function 
c'i by 

W-Jt + l \^ Pi- n -I 

and p= k+2 , , , . , N 



argmm 



2 W-Jt + 1 



p=Jt + 2 ^ 
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Given an index pair {l^f ^k^i) f ^Jk^zf - - * f Dn) need not be unique, 
resulting in the potential generation of several subgradients 
(4 . 17) . Then, 

Minimize ^^-^c-ki ( . . . ,u^) 

- £ r ch -t tuf (4,10) 

Subject To 22 ^1 i =lf • • • ijt/ 

<i,.,e{0,l} for all Vi;,,!- 



As an aside, two observations are in order. The first is that 
the search procedure needed for the computation of the relaxed 
cost coefficients in (4.9) is the most computationally 
intensive part of the entire relaxation algorithm. The second 
is that a feasible solution z^"^^-^ of a sparse problem (4.7) 
yields a feasible solution of (4.10) via the construction 

2 J 1, if z^'f'^y ...i =1 for some [ij,^^, - , i J , 
[0, otherwise. 

thus, there are generally solutions other than the one nonzero 



index solution. 
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The following Theorem 4 . 1 gives a method for evaluating 
and states that an optimal solution of (4.8) can be computed 
from that of (4.10) . Furthermore, if the solution of either of 
these two problems is G-optimal, then so is the other. The 
5 converse is contained in Theorem 4.2. 

Theorem 4.1. Let jbe a feasible solution to the two- 
dimensional assignment problem (4.10) and define w^"^^^ by 



and (i^,ij^^^)^ (0,0: 



W-Jt+l T 



and * (0,0; 

-! — •■ " p=Jt+2 " 

U W-Jt+l p. N-k+1 , P ^ n 

" ^ooi -i =0 2f Cooi ...i 2^ Uf > 0. 



(4.11) 



Then w^"^*-^ is a feasible solution of the Lagrangian relaxed 
problem and 

10 0^.;,,,(w^-^^^•u^*^...,u^) = (p.j,.,{w,; u^^^ . . . ,u^) . 

If, in addition, is optimal for the two-dimensional problem, 
then w^°^^^ is an optimal solution of the relaxed problem and 

^N-j..i ( u^*^ . . . , u^) = s_^^, ( u^^^ . . , , u^) . 
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Proof. Let w^'^""^ and be as in the hypotheses, and let 

. • ./U^) and • 

denote the objective function values of (4.8) and (4,10), 
respectively. Direct verification shows that w^"^""^ satisfies 
the constraints in (4.8) and 

0^^.^.! ( ; ; . . . , u^) = ( ; u^'' , ... , u^) . 
For the remainder of the proof, assume that is optimal for 
(4.10). Let ^'^^^ satisfy the constraints in (4.8) and define 

^li =Ei i ^rr'i -i ^or (l.,i. J ^(0,0), 
Xoo-0 if (lj^,ij^,i) = (0,0) and 

Note that satisfies the constraints in (4.10) and 



168 



CSUR.01USRI 

for any feasible solution x""**^ of the constraints (4.8) . This 
implies w""^*-^ is an optimal solution of (4.8) . Next, 

<f^-^.Au>^'', . . . , u") = (u*^". . . . , u") 

follows immediately from 

for an optimal solution of (4.10) . With the exception of one 
equality being converted to an inequality, the following 
theorem is a converse of Theorem 4.1. 

Theorem 4.2. Let u^'^"^ Jbe a feasible solution to problem (4.8) 
and define by 

^1 i = E ^fT\ i Tor (l,,i^,J^ (0,0) , 
w^^^l if (i^, i^^^) = (0, 0) and 

^N'k+l , ^ P A jr / ' • \ (4.13) 

^ooi,,,-i^ + 2^ ^0 for some (ij,,2'-' V ' 
w^Q = 0 if (Ij^,i^^^) = (0,0) and 

Then a 
feasible 

solution of the problem (4.10) and 

(pN-k.i ( ; u^^' ru^) ^ <P^.,.i ( ^ ; u^'' , . . . , U^) . 

If, in addition, i/^'^"^ is optimal for (4.8), then vj^ is an 
optimal solution of (4.10), 
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Proof: Let v/^"^"*"^ and be as in the hypotheses, and let 
0N-k^i(^~^^''/^^^^f • . . / u^) and 0V;c^i (u^;u^^^ . . . , u^) denote the 
objective function values of (4.8) and (4.10), respectively. 
Direct verification shows that satisfies the constraints in 
(4.10) and 

<pN-k.i ( w^-^^^•u^^^ . . . , u^) > 4.;c.i (w^/u^^^ . . . , u^) . 

For the remainder of the proof, assume that w^"^^^ is optimal 
for (4.8) and construct as above. Using Theorem 4.1, 
construct w^~^^^ by replacing w^'^^^ in (4.11) with w^'^^^ . Then, 
from that theorem 

( ''''''^^ . . . , u^) =0.-k.A^;u^'', . . . , u^) 

Optimal ity of w^"^""^ then implies the last inequality is in fact 
an equality. This proves the theorem. 

The Nonsmooth Optimization Problem. Next, we address the 
nonsmooth properties of the function ^^j^.j^.^^ ( u^^^ , . , . , u^) as 
explained in the following theorem. 

Theorem 4.3. (G.L, Nemhauser and L,A. Wolsey. Integer and 
Combinatorial Optimization, Section II. 3. 6. Wiley, New York, 
NY, 1988). Let 0^.^,^:^ he as defined in (4.7), let V^./c^i Cz^'^^^J 
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be the objective function value of the (N-k+1) -dimensional 
assignment problem in equation (4.7) corresponding to a 
feasible solution z^'^*^ of the constraints in (4.7), and let 
fN-k^i optimal solution of (4.7). Then, (p^_j^^^ (u^^^ , . . . , u^ ) 

5 is piecewise affine, concave and continuous in {u^"^, , . . ,u^} and 

Most of the algorithms for nonsmooth optimization are based on 
10 generalized gradients, called subgradients , given by the 
following definition for a concave function. 

Definition 4.1, At 0 = [u^^^ , . . . , u the set [u) is called 

a subdifferential of ^^.j^^i is defined for the concave 

15 function 0^_j^^^{u) by 

d0^_^,,(u) ={ge m""-^^' x...xR^-^'|^^_^^^(p^) ~0^_,,,(u) ^ (4.15) 
g^iw'U) for all we r""'*'"^ x ... x l""^"^ } , 

where go* = t\^o'= Uo'= 0 are all permanently fixed. (Recall that 
20 these were used for notational conveience only.) A vector 
g E d0^_j^^^ is called a subgradient. 
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There is a large literature on such problems, e.g., 
(J.-B. Hiriart-Urruty and C. Lemarechal . Concex Analysis and 
Minimization Algorithms I& II. Springer- Verlag, Berlin, 1993; 
K.C. Kiwiel. Methods of descent for nondif f erentiable 
5 optimization. In Lecture Notes in Mathematics 1133, A. Dold 
and B. Eckmann, eds. Springer-Verlag , Berlin, 1985; 
C. Lemarechal and R. Mifflin, Nonsmooth Optimization, Pergamon 
Press, Oxford, UK, 1978; B.T. Polyak. Subgradient method: 
A survcy of Soviet research. In C. Lemarechal and R, Mifflin, 
,==1 10 editors, Nonsmooth Optimization, pages 5-29, N.Y., 1978. 

a:: 
: ) 

Pergamon Press.; H.Schramm and J. Zowe . A version of the 
h bundle idea for minimizing a nonsmooth function: Conceptual 

i I 

idea, convergence analysis, numerical results. SI AM Journal on 

Mr, 

n Optimization, 2, No. 1 : 121-152 , February, 1992; N.Z. Shor. 

'IS 

fi 15 Minimization Methods for Non-Diff erentiable Functions . 

Springer- Verlag, Nem York, 1985; P.Wolfe. A method of 
conjugate subgradients for minimizing nondif f erentiable 
functions. Mathematical Programming Study, 3:145-173, 1975), 
and we have tried a variety of methods including subgradient 
2 0 methods (N.Z. Shor. Minimization Methods for Non- 
Diff erentiable Functions . Springer- Verlag, New York, 1985) 
and bundle methods (J.-B. Hiriart-Urruty and C. Lemarechal. 
Convex Analysis and Minimization Algorithms I& II, Springer- 
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Verlag, Berlin, 1993; K.C. Kiwiel . Methods of descent for 
nondif f erentiable optimization. In Lecture Notes in 
Mathematics 1133, A. Bold and B.Eckmann, eds. Springer -Verlag , 
Berlin, 1985) . Of these, we have determined that for a fixed 
5 number of nonsmooth iterations (e.g., twenty), the bundle 
trust method of Schramm and Zowe (H.Schramm and J. Zowe . A 
version of the bundle idea for minimizing a nonsmooth 
function: Conceptual idea, convergence analysis, numerical 
results- SIAM Journal on Optimization, 2, No. 1:121-152, 
10 February, 1992) provides excellent quality solutions with the 
fewest number of function and subgradient evaluations, and is 
therefore our currently recommended method . 

An Algorithm for Computing ^^_j^^^and a Subgradient • 
Most current software for maximizing the concave function 
15 ^u-jc+x requires the value of the function and a subgradient at 
a point ( u^""^, . . . , u ^) . Based on the previous two sections, 
this can be summarized as follows. 

1. Starting with problem (4.7), form the relaxed 
problem (4.8) ; 

20 2 . To solve (4.8) optimally, defined the two- 

dimensional assignment problem (4.10) via the transformation 
(4.9) ; 

3. Solve the two-dimensional problem (4.10) optimally; 
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4. Reconstruct the optimal solution, say ^^'^^'^ , of 
(4.8) via equation (4.9) as in Theorem 4.1; 

5. Compute the function 



5 <^*w-j..i("'^''' ...,u*') 



p=;c+2 ip=o 



■^k-^k+l -^N 



(4.16) 



10 



6. A subgradient is given by 



91 



p _ 



dul 



for i^=l, 



, and p = k+2, 



(4.17) 



Several remarks are in order. First, the optimal solution of 

the minimization problem (4.8) is required before one can 

remove the minimum sign, replace z^"^"*'-'- by the minimi zer ^'^^■^ 

and differentiate with respect to uf to obtain a subgradient 

p 

15 as in (4.17). If ij^'^"^ is an approximate solution of (4.8), 
then the subgradient and function values are only approximate 
with accuracy depending on that of (J^"^""^ . Although one can 
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evaluate the sums in (4.16) and (4.17) in a straight forward 
manner, another method is based on the following observation. 
Given the multiplier vector • • ' -'-^^ 

{(■^jt^-^/c.i) ' ^jt.i(^^.i) )}^'' enumeration of indices jl^. ik.i} of 

5 (or the first 2 indices of w^'^^'^ constructed in equation 

(4.9)) such that for (1^(1,,,) , i^,, (I,,,) ) ^ (0, 0) 

and (1^(1^^^), ij^^^d^^^)) = (0,0) for I^^^ = 0 regardless of whether 
p/o^Q = l or not. (The latter can only improve the recovered 
feasible solution.) The evaluation of the bracketed quantity 
10 in (4.17) for a specific 1 and p = k-\-2,...,N is one minus 
the number of times the index value ip appears in the (p-k+l)^^ 
position of the ("N-Jc+l J - tuple in the list 

Finally, we have presented a method for computing one 
15 subgradient . If i)^'^^-^ is the unique optimal solution of (4.8) , 
then is dif f erentiable at u, g is just the gradient at 

u, and the subdif f erential 8<?^_^^^ ( u) = { g} is a singleton. If, 
corresponding to u, 0^'^^^ is not unique, then there are 
finitely many such solutions, say { i^^'^^^l) / - - - / (m) } with 

20 respective subgradients {g(l) , ... ,g(m) }, the subdif f erential 
9<^(u) is the convex hull of (g (1) , . . . , g (m) } (J.L. Goffin. On 
convergence rates of subgradient optimization methods, 
mathematical programming. Mathematical Programming, 13:329- 
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347, 1977). These nonunique solutions arise in two ways. 
First, if the optimal solution of the two-dimensional 
assignment problem is not unique, then one can generate all 
optimal solutions of the two-dimensional assignment problem 
5 (4.10). Corresponding to each of these solutions and to each 
index pair [l^,i^^^)±r). each solution, compute the indices 
( j^^2' • • • '^w^ in (4.9). If these j's are not unique, then we 
can enumerate all the possible optimal solutions v^'^^'^ of 
(4.8). Given these solutions, one can then compute the 
10 corresponding subgradients from (4.17). In most nonsmooth 
optimization algorithms, one uses an G-subdif f erential . 

Definition 4.2, At u = (u^^^ , . . . ,u^) the set d^0^_j^^^(u) is 
called an e-subdiff erential of ^^.^^^ ^^d is defined for the 
15 concave function by 

g^{w-u) + e for all w 6l"*-^*^x . . . x r""*^}, 

where qQ=WQ=UQ=0 are all permanently fixed. (Recall that 
these were used for notational convenicence only. ) A vector 
g 6 d^(P^_i^^j^(u) is called an e-subgradient. 
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IV. 3. 3 The Recovery Procedure. 
The next objective is to explain a recovery procedure, i.e., 
given a feasible (optimal or suboptimal) solution of (4.10) 

(or w^'^""^ of (4.8) constructed via Theorem 4.1), generate a 
feasible solution z^'^^^ of equation (4.7) which is close to 
in a sense to be specified. We first assume that no variables 
in (4.7) are preas signed to zero; this assumption will removed 
shortly. The difficulty with the solution v^'^^^ in Theorem 4.1 
is that it need not satisfy the last (N-k-1) sets of 
constraints in (4,7) . (Note, however, that if is an optimal 
solution (4.10) and w^'^""^, as constructed in Theorem 4.1, 
satisfies the relaxed constraints, then w^'''^'"^ is optimal for 

(4.7) .) The recovery procedure described here is designed to 
preserve the 0-1 character of the solution of (4.10) as far 
as possible: If wf . =1 and 1,^0 or i.^.^0, the 

corresponding feasible solution z^'^^^ of (4.7) is construed so 
that = ^ for some (ij^^2f - - - / ^n) - this reasoning, 

variables of the form Znn^"^^ can be assigned a value of one 

2 

in the recovery problem only if Wqq = 1 , However, variables 
Zooif^l-i^ will be treated differently in the recovery procedure 
in that they can be assigned either zero or one independent of 
the value w^q . This increases the feasible set of the 
recovery problem, leading to a potentially better solution. 
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Let {{l^{l^^^),i^^^{l^^^)Y^"'^_ be an enumeration of indices 

{1^, of ^ (or the first 2 indices of v/^'^^^ constructed in 

equation (4.11)) such that 

^1 ii 1 i (I 1 ~ 1 for 

Vi(^..i)) = {0'0) for I,.,=0 
regardless of whether w^^^l or not. To define the (N-k) - 
dimensional assignment problem that resores feasibility, let 

= n c^"^ = for k=l 

^ Ci J ^ — Ci fiwMw i (4.18) 

: = C • • ■ ' * • ' 

ffl ^1 f ^2. . .;t+i^ -^2 ^ ^2. . -^3 ^-^a. . • • * -^Jt ^ -^JtJt+l' -^it+i ^ -^Jt+i ' -^Jt+2 • • • -^w 

Tor A: ^ 2 and I . = 0 , . . . , L. , 



rii 

in 



where 



10 for m = 2, , . and where "o" denotes function composition. 
For example, I^j = ^2(^3) 1234= 12^13(14)== 12(13(14)). 

Then the (2\r-Jc>) -dimensional assignment problem that 
restores feasibility is 
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, N-k N-k 

Minimize J, Z; 

Subject To E . ^f;^^^^ . . . = 1 ' -^fc.i = 1 ' • • • ' -^it.i ' 



(4.20) 



for 2^= 1, . . . and p = A-+3, . . . , 



E^"^ - 1 ' = 1 



<;^i,.,...i,^{0,l} for all 

Let Y be an optimal or feasible solution to this {N-k) - 

dimensional assignment problem. The recovered feasible 
solution is defined by 



1, if i, = i,(l2...;,,i), i2 = i2(^2...;:.l)' 

i3 = i3(^2...;c.i)'--- "^k-^k^^kk^i"^ ' (4,21) 
0, otherwise. 



Said in a different way, the recovered feasible solution 
corresponding to the multiplier set { u^"^, . . . , u^} is then 
defined by 
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1, if , =1, 

0, otherwisSf 



where Im^-k+i is defined in (4.19) and "o" denotes function 
composition . 

The next objective is to remove the assumption that all 
5 cost coefficients are defined and all zero-one variables are 
free to be assigned. We first note that the above 
construction of a reduced dimension assignment problem (4.20) 
is valid as long as all cost coefficients c^''^ are defined and 
all zero-one variables in z^'^ are free to be assigned. 

10 Modifications are necessary for sparse problems. If the cost 
coefficient c^'n'^.^ n . is well-defined and the zero-one 

variable zTn^ m n m i is free to be assigned to zero or 
one, then define cj^ = c^'n'^ m n m as in (4.18) with 

Zj^'^. being free to be assigned. Otherwise, zf is 

15 preassigned to zero or the corresponding arc is not allowed in 
the feasible set of arcs. To ensure that a feasible solution 
exists, we now only need ensure that the variables z^^'^^o-o 
free to be assigned for l^^j^ =0 , 1 , ... , L^^^ with the corresponding 
cost coefficients being well-defined. (Recall that all 
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variables of the form zj^olo^ and ^q'q^^q.^q are free (to be 
assigned) and all coefficients of the form c^^^^q and 
^o'oipO-o well-defined for Ij^ =0,..,,Lj^ and ip =0,.../Mp for 

p=k'^l, ...,N.) This is accomplished as follows: If the cost 
5 coefficient cf;(V^;^ji^^^(,^^^jo.-.o is well-defined and zf;(tj^,,^^^(,^^^jo.-.o 
is free, then define cf^^J'o-o = ^UCo^.^.ih.i^o.-^o ^ith Zi'^^^o-o being 
free. Otherwise, since all variables of the form z^^Q^^^Qand 
^oi^lo'-o free to be assigned with the corresponding costs 

being well-defined, set cf^"^fo.-.o = cf,u,,,)oo--.o + Cot,^^fv,)o-.-o / where the 
10 first term is omitted if Ij^(lj^^j^) =0 and the second, if 
^k+idk+i)^^ ' For l]^(lk+i)^0 and ij^^idk+i) =0 , define c^'q = c^'^"^ . 



IV. 3. 4. The last step k = N-2 



0,^Minimize 2^ 2^ c; ^ ^^VAz^) 



2 



2 



subject To 1. <.,i„=l' lw-i = l'-' Vi' (4.22) 

e {0,1} for all l^.,,i^. 

The description of the algorithm ends with k = The 

15 resulting recovery problem defined in (4.10) with k = N-2 is 
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Let Y be an optimal or feasible solution to this two- 
dimensional assignment problem. The recovered feasible 
solution is defined by 



or if ( Vi.iJ = (0,0) , 
0, otherwise. 



Said in a different way, the recovered feasible solution is 
then defined by 



1, if y, . =1, (4.24) 

0, otherwise. 



where Im-N-i is defined in (4.19) and "o" denotes function 
composition. To complete the description, let 

, i^{l^)V\ be an enumeration of indices of Y 

such that i',^.^„^,,,^„^,=l for , ) MO, 0) and 

(l^.^ (Ij,) , ij^(Ij^))= (0, 0) for 1^=0 regardless of whether 3^00=1 or 
not. Then the recovered feasible solution can be written as 



9". . . =< 



0, otherii^ise. 



(4 .25) 



IV. 3. 5. The Upper and Lower Bounds 
The upper bound on the feasible solution is given by 
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and the lower by 0^(u'^ , ..,,u^) for any multiplier value {u^,..,,u^) . 
In particular, we have 

cP^(u^...,u^) ^ ^^(u^...,u^) ^ ^^.(i"'')^ ^^,(^''). (4.27) 

where {u^,.,.,u^) is any multiplier value, (u^,.../U^) is a 
maximizer of ^^(u^, u^) , is an optimal solution of the 

multidimensional assignment problem (4.7) and is the 

recovered feasible solution defined by (4.25). Finally, we 
conclude with the observation that V^(z)=V2{Y) where Y is the 
optimal solution of (4.23) used in the construction of z in 
equations (4 . 24) - (4 .26) . 

IV. 3. 6. Reuse of Multipliers 
Since the most computationally expensive part of the 
algorithm occurs in the maximization of ^^.jc+it the development 
of algorithms that can make use of hot starts or multipliers 
close to the optimal are fundamentally important for real-time 
speed. The purpose of this section is to demonstrate that the 
multiplier set obtained at stage k^l provide good starting 
values for those obtained at step k-hl . 
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Theorem 4.4, Let (u^,...,u^) denote a multiplier set obtained in 
the maximization of 0j^{u'^ , ,..,u^) . Then this multiplier set 
satisfies the string of inequalities 

0^{u\ . . .,u^) ^ V,{u\ . . .,u^) ^ . . . ^^^{u^) ^0^^V^{J) , (4.28) 

where after the first step in the maximization of 0^ the 
5 multipliers are not changed in the remaining steps. 
Furthermore, to improve this inequiality, let 
% (u2^^'^-^^S...,u^'^-^^^) denote a maximizer of ( u^^^^ u ^) . 

Then, we have 

^ <5„_j,^, (u2^*'"-^^S . . . , u'^-^-'^'^) (4.29) 



Li 



Proof: Suppose we have a value of the multipliers {u^'^^ , ... , Uj^) 
10 obtained in the maximization of 

^i^^^^^e ^„.j^,,(z*'-'^^S-u**2^ . . .,u") (4.30) 
Subject To Yl ZiT\ i =1, l.= l, . . .L., 



E N-'cn . =1 i =1 M 



where 
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EW-Jt+1 W-Jt + 1 

^1 i i ^1 i i 



p=k+2 ip=0 



(4.31) 



w-Jt+i ^ p w-Jt+i _ ^ Y*^ p 



p=Jt+2 



p=;:+2 ip=0 



These need not be the maximizer; however, we do assume that we 
have solved the above minimization problem optimally to 
evaluate (u^^^, u ^) . Just as in the definition of the 

recovery problem discussed earlier, let 
5 , ij^^^ (1^^^) )}2^^|=o be an enumeration of indices {Ijcrik+i) 

of the optimal solution of problem (4.10) (or the first 2 
indices of the solution u^"^"**"^ of the relaxed problem (4.8) 
constructed in equation (4.11)) such that ^ . . .=1 for 

■^k^-^k^l' * ^k+l^-^k^V 

U;c(Vi)'i..i(Vi))'' (O'O) and V,) , i,,, (1,,,))= (0, 0) for 1^,^=0 
10 regardless of whether ^l^-l or not. 

Then, the final value of the objective function can be 
expressed as 
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Minimize ( z"-**^- u'^^^^ u ^) (4.32) 

Subject to Y.. <u;;,)i,,,(i,.,)Wi„=l'-^^.i = l---^..i' 



where 



EW-Jt+1 ^ f 

£ £ < 



If now another constrain set, say the (A:+2)^^ set, is added to 
5 this problem, one has 

^^=jt.i ^ Minimize ^^.^^^^ (z^-^^^• u^^^^ u ^) ^ (4.34) 

Minimize E . <"(t;,)i,.,(2,.,)...i/ £ < 

- £ £ "f. 

p=J:+2 ip=0 " 
^"'^^'^^^ to E. <'(t;,)i..,(i..,)W.i„=l' = 

- 1 - 1 M 



186 



CSUR.01USRI 

Since the constraint /J n'^ m a m ..i -1 fo^ 

ij^^2~'^ t - is imposed, the feasible region is smaller, so 

that one has 

where the fact that ^^^j^^j^ does not depend on the multiplier set 
u^"*"^ due the added constraint set. Also, the last equivalence 
follows from the fact that ^i^^;c+i ( ^^^^ / - / is the relaxed 
10 problem (with k replaced by k-\-l) for the recovery problem in 
step k of the above algorithm. Continuing this way, one can 
compute (u'^,...,u^) at the first step {k=l) , fix them thereafter, 
and perform no optimization at the subsequent steps to obtain 

15 

where 02 is defined in (4.22). 

To explain how to improve this inequality, let 
(u^^^'^^-^^S^.^u^'^-^^^) denote a maximizer of <^^_^^^ ( u^^^^ u ^) . 
2 0 Then by the same reasoning one has the result as stated in the 
theorem. 
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V. SUMMARY OF THE LAGRANGIAN RELAXATION PROBLEM 
Given the multidimensional assignment problem 



Minimize \\ c, . z. . 

Subject To Y\ z, , = 1 (i = 1, . . .M, ) , (5.1) 
V3...i, 

.E (ip = l, . . .Wp and p = 2, . . . , W-1) , 



•-1^2 • • '-^N-l 



z. .6(0,1} for all 



where C(J^.o is arbitrarily defined to be zero and is included 
for notational convenience and where the superscript N on both 
5 c and z is not an exponent, but denotes the dimension of the 
subscripts and the assignment problem as stated in (3.5). In 
relaxed and recovery problems Co^.o need not be zero! In this 
problem, we assume that all zero-one variables ^i^-i^ with 
precisely one nonzero index are free to be assigned and that 
10 the corresponding cost coefficients are well-defined. (This 
is a valid assumption in the tracking environment (A.B. Poore . 
Multidimensional assignments and multitarget tracking, 
partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
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Computer Science, American Mathematical Society, Providence, 
R.I., 19:169-198, 1995; A.B. Poore . Multidimensional 
assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
5 Computational Optimization and Applications, 3:27-57, 1994).) 
Although not required, these cost coefficients with exactly 
one nonzero index can be translated to zero by cost shifting 
(A.B. Poore and N. Rijavec. A lagrangian relaxation algorithm 
for . multidimensional assignment problems arising from 

10 multitarget tracking. SIAM Journal of Optimization, 3, No. 
3:544-563, 1993) without changing the optimal assignment. 

Having explained many of the relaxation features, it is 
now appropriate to present the Lagrangian relaxation 
algorithm, which iteratively constructs a feasible solution to 

15 the N-dimensional assignment problem (5.1). 

Algorithm 5.1 (Lagrangian Relaxation Algorithm) To 
construct a high quality feasible solution, denoted by z^, of 
the assignment problem (5.1), proceed as follows: 

0. Initialize the multipliers (u^*^, u ^) , e.g., 
20 (u^*2^...,u^) = (0,...,0) . 

For k-l,,„,N-2, do 

1. For the Lagrangian relaxed problem (4.8) from the problem 
(4.7) by relaxing the last {N-k-1) sets of constraints . 
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2, Use a nonsmooth optimization technique to solve 

Maximize { ^^.^^^ J u^eR^^'^ for p = k^2, ...,N (5.2) 

with Uq=0 being fixed], 

where 0^_^^^[u^^'^ , ... , u^) is defined by equation (4,8), The 
algorithm in Section IV, 3, 2 provides one method for 
computing a function value and a subgradient out of the 
5 subdifferential at ( u^^^, u ^) , as required in most 

nonsmooth optimization techniques, 
3 . Given an approxiimte orcptiiml maxiirdzer of (5,2) , say [Q^^^ , - ,(i^), 
let denote the optimal solution of the two-dimensional 
assignment problem (4,10) corresponding to this maximizer 
10 of 0^_^^^{u^'\...,u'') , 

4, Formulate the recovery {N-k) -dimensional problem (4.19), 
modified as discussed at the end of subsection IV, 3, 3 for 
sparse problems . At this stage, z^ , as defined in (4,21), 
contains the alignment of the indices {ii/ ik+i) - 
15 Formulate the final two-dimensional problem (4.22) and 
complete the final recovered solution as in (4.25). 

To explain the lower and upper bounds, let 0^ be as 
defined in (4.8) with k=l, let V^{z^) be the objective function 
value of the 2^- dimensional assignment problem in equation 
20 (3.5) corresponding to a feasible solution of the 
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constraints in (3.5), and let ^'^'^'^ be an optimal solution of 
(3,5) . Then, 

^^(u^...,u^) ^ \/^(z^) ^ V^i^"") (5.3) 

is the desired inequality. 
COMMENTS : 

1. Step 2 is the computational part of the algorithm. 
Evaluating ^i^.jc+i and computing a subgradient used in the 
search procedure requires 99% of the computing time in 
the algorithm. This part uses two-dimensional assignment 
algorithms, a search over a large number of indices, and 
a nonsmooth optimization algorithm. It is the second 
part (the search) that consumes 99% of the computational 
time and this is almost entirely parallelizable . 

2. In track maintenance, the warm starts for the initial 
multipliers for the first step are available. For the 
relaxation procedure, initial multipliers are available 
in step 2 from the prior step in the algorithm. 

3. There are many variations on the above algorithm. For 
example, one can compute a good solution on the first 
pass {k=l) and not perform the optimization in (2) 
thereafter. This yields a great solution. Thus one can 
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continue the optimization at the first pass, and 
immediately compute quality feasible solutions to the 
problem. 

VI. LAGRANGIAN RELAXATION ALGORITHM 
FOR THE THREE-DIMENSIONAL ALGORITHM 

In this section, we illustrate the relaxation algorithm 

for the three-dimensional assignment problem. Having 

discussed the general relaxation scheme, 



M, 



Min 



imize £ £ £ 

1 =n 7 =n 7 =n 123 123 



ii=0 i2=0 ^3=0 



Subject To 



1^ = 0 13=0 ^'^^ ' 



(6.1) 



il=0 i3=0 12-3 



1^=0 12=0 ^ 2 3 



z . . . 6(0,1) for all 1,12^2' 



To ensure that a feasible solution of (6.1) always exists for 
a sparse problem, all variables with exactly one nonzero index 
(i.e., variables of the form ^i^oo'^a^o^ ^ooi^ ^^-^ ip=l,...,Mp 

and p-1,2,3 are assumed free to be assigned with the 
corresponding cost coefficients being well-defined) . This 
assumption is valid in the tracking environment (A.B. Poore. 
Multidimensional assignments and multitarget tracking, 
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partitioning data sets, I.J. Cox, P. Hansen, and B. Julesz, 
editors, DIMACS Series in Discrete Mathematics and Theoretical 
Computer Science, American Mathematical Society, Providence, 
R.I,, 19:169-198, 1995; A.B. Poore . Multidimensional 
5 assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994). 

VI . 1 . The Laqrangian Relaxed Assignment Problem 
10 The three-dimensional assignment problem (6,1) has three 

sets of constraints and one can describe the relaxation by 
relaxing any of the three sets of constraints, the description 
here is based on relaxing the last set of constraints. A 
(M3+I) -dimensional multiplier vector (i.e., ^3 ^ M^^^"^^ 
15 associated with the 3"^^ constraint set will be denoted by 
=\uq , u^ , ... , u^^^'^ with Uq=0 being fixed for notational 
convenience. (The zero multiplier Uq=0 is used for notational 
convenience.) Now, the three-dimensional assignment problem 
(6.1) is relaxed to a two-dimensional assignment problem by 
20 incorporating last set of constraints into the objective 
function via the Lagrangian. Although any constraint set can 
be relaxed, we choose the last set of constraints for 
convenience. The relaxed problem is 



193 



CSUR.01USRI 




(6.2) 



*^(u^) = Minimize (^^(z^/u^) 



= Minimize 



Subject to 



One of the major steps in the algorithm is the 



maximization of ^^{u^) with respect to the multiplier vector 
. It turns out that 0^ is a concave, continuous, and 
piecewise affine function of the multiplier vector u^, so that 
the maximization of ^3 is a problem of nonsmooth optimization. 
Since many of these algorithms require a function value and a 
subgradient of ^3 at any required multiplier value (u-^) , we 
address this problem in the next subsection. We note, 
however, that there are other ways to maximize 0^ and the next 
subsection addresses but one such method. 

VI . 2 . Properties Lagrangian Relaxed Assignment Problem 
For a function evaluation of ^3, we show that an optimal 
(or suboptimal) solution of this relaxed problem (6.2) can be 
constructed from that of a two-dimensional assignment problem, 
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Then, the nonsmooth characteristics of are addressed, 

followed by a method for computing a subgradient. 

Evaluation of ^3. Define for each [i^, i^) an index 

73 = 73(2.^,12) and a new cost function c?^^^ by: 

J 3 3 , I (^-3) 

73^73(2^,12) =arg/nin|Ci^i^i^ + Ui^ | 23 = 0, 1, ...,M3|, 



2 












2 


-t 








i3=0 



Minimum 



Given an index pair (2^,22) , 73 need not be unique, 
resulting in the potential generation of several subgradients 
(6-11), (We further discuss this issue at the end of the 
Subsection. ) 

Now, define 

3 



0^{u') -Minimize 0^{z^;u')^2l 

1^=0 12=0 ^ ^ ^ ^ i3=0 ^ 

Subject To 22 ^li=^r 2^ = 1,... ,M^, (6.4) 

1^ = 0 

2 _ . _ 

2-^ ^i^ij ^' ^2 1/ • • • '^2' 

. 6 { 0, 1 } for all i.fi^. 

As an aside, two observations are in order. The first is that 
the search procedure needed for the computation of the relaxed 
cost coefficients in (6.3) is the most computationally 
intensive part of the entire relaxation algorithm. The second 
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is that a feasible solution of a sparse problem (6.1) yields 
a feasible solution of (6.4) via the construction 

lo, otherwise. 

5 Thus the two-dimensional assignment problem (6.4) has feasible 
solutions other than those with exactly one nonzero index if 
the original problem (6.1) does. 

The following Theorem 6.1 states that an optimal solution 
of (6.2) can be computed from that of (6.4) . Furthermore, if 
10 the solution of either of these two problems is e-optimal, 
then so is the other. The converse is contained in Theorem 
6.2 . 

Theorem 6.1. Let Jbe a feasible solution to the two- 

m dimensional assignment problem (6.4) and define by 

15 3 2... (^-5) 

^3 = h (ii.i2) ^ (0,0) , 

^V2i3 = 0 ^3^J3 (ii/ij) ^ (0/0) / 

^ooi3 =1 if Cooi3 + "i3^0, 
^ooi3 =0 if Cooi3^"i3>0. 

Then is a feasibl e sol u tion of the Lagrangi an relaxed 
problem (6.2) and 0^{w^ ; u'^) = 0^{w^ ; u-^) . If, in addition, is 



m 
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optimal for the two-dimensional problem, then is an optimal 
solution of the relaxed problem and ^^(u"^) =<^(u^) . 

With the exception of one equality being converted to 
5 an inequality, the following theorem is a converse of 
Theorem 6.1. 

Theorem 6 •2. Let jbe a feasible solution to problem (6.2) 
and define by 

^^L=S for {i^,i^} (0,0) , 

23=0 

w^Q = 0 if (i^, i^) = {0,0) and Cqq^^ + uI^>0 for all i^, 
Wqq = 1 if (i^, i^) = (0, 0) and Cqq^^^-^ u^^<, 0 for some i^. 

Then is a feasible solution of the problem (6.4) and 
0^{w^ } u-^) ^ (p^[w^ ; u-^) . If, in addition, is optimal for (6. 2), 
then is an optimal solution of (6.4), 0^{w-^ ; u^) = <p^{w^ ; u^) 
and 0^{u^) =0^(u^) . 

An Algorithm for Computing 0^ and a Subgradient. 

Given the problem of computing (P^iu^) and a subgradient 
of ^^(u^) can be summarized as follows. 



10 



15 
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1. Starting with problem (6.1), form the relaxed problem 
(6.2) . 

2. To solve (6.2) optimally, define the two-dimensional 
assignment problem (6.4) via the transformation (6.3). 

3. Solve the two-dimensional problem (6.4) optimally. 

4. Reconstruct the optimal solution, say , of (6.2) via 
equation (6.5) as in Theorem 6.1. 

5 . Then 



M, W, 



1^=0 i2=0 i3=0 ^ ^ ^ ^ ^ ^ 

£ ^^i, £ £ ^luu-^ 



23=0 L 1^-0 22=0 



(6.7) 



6. A subgradient is given by substituting v^^ into the 
objective function (6.2), erasing the minimum sign, and 
taking the partial with respect to u^^ . The result is 



2,=0 2, = 0 ^ ^ ^ 



for i^ = l. 



(6.8) 



VI . 3 . The Recovery Procedure 
The next objective is to explain a recovery procedure , 
i.e., given a feasible (optimal or suboptimal) solution of 
(6.4) (or of (6.2) constructed in Theorem 6.1), 

generate a feasible solution of equation (6.1) which is 
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close to in a sense to be specified. We first assume 
that no variables in (6. 1) are preassigned to zero; this 
assumption will be removed shortly. The difficulty with the 
solution vr" in Theorem 6.1 is that it need not satisfy the 
5 third set of constraints in (6.1). (Note, however, that if 
is an optimal solution for (6.4) and w^, as constructed in 
Theorem 6.1, satisfies the relaxed constraints, then is 
optimal for (6.1).) The recovery procedure described here is 
designed to preserve the 0-1 character of the solution of 
10 (6.4) as far as possible: If ^i^i^^l and i^^O or ia^O, the 
corresponding feasible solution of (6.1) is constructed so 
that ^i^i^i^ = 1 fo^ some (11,22,23). By this reasoning, variables 

of the form Zg^oij *^a^ assigned a value of one in the 

2 3 
recovery problem only if . However, variables ^0013 

15 be treated differently in the recovery procedure in that 
they can be assigned either zero or one independent of the 
value . This increases the feasible set of the recovery 
problem, leading to a potentially better solution. 

Let , i^[l^)Y'^ ^ be an enumeration of indices {11,23} 

20 of (or the first 2 indices of constructed in equation 
(6.5) ) such that =1 for (i^ds) ^^2(^2))'' 

, [Q ,0) for l2=0 regardless of whether w^o^^^ 
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not. (The latter can only improve the quality of the feasible 
solution. ) 

Next to define the two-dimensional assignment problem 
that restores feasibility, let 

2 _ 3 _ (6.9) 

Then the two-dimensional assignment problem that restores 
feasibility is 



Minimize ^ £ ^tu^lu 



(6.io; 



2,=0 i,=o " 



Subject to 22 "1' I2 = l/...fi2/ 



6(0,1} for all l^,i 



2/ -3' 



The next objective is to remove the assumption that all 
cost coefficients are defined and all zero-one variables are 
free to be assigned. We first note that the above 
construction of a reduced two-dimensional assignment problem 
(6.11) is valid as long as all cost coefficients are defined 
and all zero-one variables in are free to be assigned. 
Modifications are necessary for sparse problems. If the cost 
coefficient (12)12 (22)^3 well-defined and the zero-one 

variable (12)^2(12)^3 ^® free to be assigned to zero or one, then 
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define cl^i^ = cl^^.^^^^^^^^^^ as in (6.9) with zl^^^^^^^^^^^^^ being free 
to be assigned. Otherwise, ^i^i^ preassigned to zero or the 
corresponding arc is not allowed in the feasible set of arcs. 
To ensure that a feasible solution exists, we now only need 
5 ensure that the variables z^^^ are free to be assigned for 
1^ = 0 ,1, ...f with the corresponding cost coefficients being 
well-defined. Recall that all variables of the form zl^^^, ^oi^o 
and Zo^oig free (to be assigned) and all coefficients of the 

corresponding form c/^qq, Cq^^^q and c^^i^ need to be defined. This 
10 is accomplished as follows: if the cost coefficient cl^^^^^^^^j^^^^^ 

is well-defined and ^1^(12)22(12)0 free, then define 

23 2 

with z_^^Q being free. Otherwise, since all 

variables of the form zl^^ and Zq^^q free to be assigned 

with the corresponding costs being well-defined, set 

2 3 3 

15 ^i^{i2)oo'^ f where the first term is omitted if 21(12) =0 

and the second, if 12(12) =0. For ii(l2)=0 and i2(l2)=0 / define 
^00 ~ ^000 • 

VI. 4, The final recovery problem 
The recovery problem for the case N = 3 is 
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(6.11) 



i,=o i,=o ^ ^ ^ ^ 



Subject to 2J =1/ Ij = 1/ •••/ 

i3 = 0 ^ ' 



1^=0 
2 



z^^i^ e {0,1} for all l^,i^. 

Let y be an optimal or feasible solution to this two- 
dimensional assignment problem. The recovered feasible 
5 solution is defined by 

i=0 



y (612) 
n ^i^2^3 [0, otherp/ise. 

n Said in another way, let { ( -^2 ^ -^3^ ' -^3 ^ -^3^ }i^=o enumeration 



Ti of indices {12,13} of Y such that i", , , . . =1 for 

2^ 3 ' ' "^3 ' 3 ' 



10 (l^dj) ,13(13) ) ^ (0,0) and (12(13) ,l3(l3)) = (0,0) for 13 = 0 
regardless of whether i^oo=l or not. Then the recovered 
feasible solution can be written as 

^3 ^/1, if i^ = i^{l^^) , i^ = i^(l^^) ,i^ = l^{l^) , ^^^^^ 
^1-^2^3 [0, otherwise. 

15 

The upper bound on the feasible solution is given by 

^ ^ ^ 

VA2')=i: E E c's' (6.14) 



ii=o 1^=0 13=0 ^ ' ^ » ^ ^ 



20 and the lower by ^^(u-') for any multiplier value (u^) 
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In particular, we have 

where is any multiplier value, is a maximizer of ^3(u^), 
—3 is an optimal solution of the multidimensional assignment 
5 problem and is the recovered feasible solution defined by 
(6.13) . 

VII. OTHER RELAXATIONS FOR THE MULTIDIMENTIONAL 
ASSIGNMENT PROBLEM 

In this section, we briefly describe other possible 
10 relaxations and their implications. These include linear 
programming relaxations and the corresponding duals. 

Recall from Section II that one starts with the 
definition of the zero-one variable z. . =1 and cost 
coefficients to form the problem 

Minimize c. . z, . 

Subject To 53 ^i...i=l {i^ = l,-,M^) , 

E ^i,...i=l {i, = l,-,M,), (7.1) 



'1 -"w 



2 0 ^rVi^i-^N 



(ip=l,-,Mp and p = 2,-,N-l) , 
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where Co.o is arbitrarily defined to be zero. Here, each group 
of sums in the constraints represents the fact that each non- 
dummy report occurs exactly once in a "track of data". One 
can modify this formulation to include multiassignments of 
5 one, some, or all the actual reports. The assignment problem 
(3,5) is changed accordingly. For example, if z/^ is to be 
assigned no more than, exactly, or no less than times, then 
the " = 1" in the constraint (3,5) is changed to " ^ - , n^^, " 
respectively. Modifications for group tracking and 

10 multiresolution features of the surveillance region will be 
addressed in future work. In making these changes, one must 
pay careful attention to the independence assumptions that 
need not be valid in many applications. 



15 and O.E, Drummond. Track initiation and maintenance using 
multidimentional assignment problems. In D.W, Hearn, W.W, 
Hager, and P.M. Pardalos, editors, Network Optimization , 
volume 450, pages 407-422, Boston, 1996. Kluwer Academic 
Publishers B.V.) significantly extends the formulation of the 

20 track maintenance and initiation to new approaches. The 
discussions of this section apply equally to those formu- 
lations . 



Again, the recent work of Poore and Drummond (A.B. Poore 



VII . 1 . The Linear Programming Relaxation 
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In the linear programming relaxation, one replaces the 
zero-one variables constraints 

z. . e{Q,l} for all (7.2) 
with the constraint 
5 O^z. . ^ 1 for all (7.3) 

Then, the problem (7.1) can be formulated as a linear 
programming problem with the constraint (7.2) in (7.1) 
replaced by (7.3) with a special block structure to which the 
Dantzing-Wolf e decomposition is applicable. Of course, after 
10 solving this problem, one must now recover the zero-one 
character of the variables in (7.1) and there are many ways to 
do this, such as using the two-dimensional assignment 
problems. Commercial software is also available. 

15 VII. 2. The Linear ProgramminQ Dual and 

Partial Lagrangian Relaxations 

Given the linear programming relaxation, one can formulate 

the dual problem or the partial Lagrangian relaxation duals 

with respect to any number of constraints. In particular, 

20 this is precisely what is done in Section 3 on the Lagrangian 

relaxation algorithm presented. The much broader class of 

algorithms provided in the US Patent 5537119 (Aubrey B. Poore, 

Jr., Method and System for Tracking Multiple Regional Objects 

by Multi -Dimensional Relaxation, US Patent Number 5537119, 
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filed 14 March 1995, issue date of 16 July 1996) can be 
modified to remove the zero-one character when one relaxes M 
sets of constraints to an {N-M) -dimensional problem and 
recovers with an {N-M +1) -dimensional problem. This avoids 
5 the problems associated with the NP-hard {N-M)- and {N-M +1) - 
dimensional problems. However, to restore the zero-one 
character, one can do it sequentially with an assignment 
problem or with one of the many zero-one rounding techniques. 
These formulations are easy to work out and thus the details 

10 are omitted. 

The complete dual problem is another way of solving the 
LP problem and may indeed be more efficient in certain 
applications. In addition, the solution of this dual problem 
may provide excellent initial approximations to the 

15 multipliers for Lagrangian relaxations. 

VIII. HOT STARTS FOR TRACK MAINTENANCE 

AND INITIATION: BUNDLE MODIFICATIONS 

Thus reuse of multipliers and the first proof that this 

20 reuse is actually valid was presented in this section. The 

reuse in track maintenance is demonstrated and discussed in 

the work of Poore and Drummond (A.B. Poore and O.E. Drummond. 

Track initiation and maintenance using multidimensional 

assignment problems. In D.W. Hearn, W.W. Hager, and P.M. 
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Pardalos, editors, Network Optimization, volume 450, pages 
407-422, Boston, 1996. Kluwer Academic Publishers B.V.). The 
only item left is the modification of the bundle of 
subgradients for the use with the multiplier values as one 
goes through the recovery problem as in Section IV. 3. 6 or as 
one moves the window in track maintenance. It is this aspect 
of the nonsmooth optimization that adds an order of magnitude 
to the speed of the algorithms. This work is based on both a 
mathematical proof as in Section IV. 3. 6 as well as 
computational experience and heuristics. 

IX. ADAPTIVE AND VARIALBE WINDOW SIZE SELECTION 
These data association algorithms , based on the 
multidimensional assignment problem, range from sequential 
processing to many frames of data processed all at once! The 
data association problem for multisensor-multitarget tracking 
is formulated using a window sliding over the frames of data 
from the sensors. This discussion focuses on the work of 
Poore and Drummond and not on the formulations in Section III. 
Firm data association decisions for existing track are made 
at, say frame k, with the most recent frame being k+M. 
Decisions after frame k are soft decisions. Reports not in 
the confirmed tracks are used to initiate tracks over frames 
numbered k-I to Jc+M. 
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The length of these windows varies from sequential 
processing, which corresponds precisely to I = 0 and M = 1, to 
user defined large values of J and M. There is a marked change 
in performance over this range. As the two windows become 
5 exceedingly long, there is little statistically significant 
information gained and indeed performance can degenerate 
slightly. Thus, one must manually find the correct window 
lengths for performance in a given scenario and the algorithms 
do not change for a changing environment. Thus we propose an 
10 adaptive method for adjusting the window lengths. (The method 
has been highly successful for selecting the correct window 
length for multistep methods in ordinary differential 
equations . ) 

1. Compute the solution for different window lengths that 
15 overlap and differ by one or two frames. 

2. Compare the solution quality, i.e., the value of the 
objective function, for two adjacent window lengths. 

20 3. If there is a predefined improvement in either direction, 
we then, for stability, repeat the exercise for a shorter 
or longer on the first try. If there is consistent 
information, we adjust the window size in the indicated 
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direction. This can be given a well defined mathematical 
formula in terms of the assignment problems of different 
dimensions . 

5 



10 

X. SENSITIVITY ANALYSIS 
For sequential processing in which the two-dimensional 
assignment algorithm is used, one can use the LP sensitivity 
analysis theory and easily obtain the corresponding answers. 

15 For the multif rame processing, the optimal solution is not 
available; however, there are still two approaches one can 
use. First, the basic algorithm can perform the same 
sensitivity analysis at each stage (loop) of the algorithm as 
is done in the two-dimensional assignment problem since the 

2 0 evaluation of function ^^.j^^^^ is equivalent to a two-dimensional 
assignment problem. Alternately, one could use an LP 
relaxation of the assignment problem and base the sensitivity 
on the resulting LP problem. We currently see this as an 
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important step in finding even better solutions to the 
assignment problem if so desired. 

5 



10 

XI . NEW AUCTION ALGORITHMS , 
In this section we present a new auction algorithm for the 
two-dimensional assignment problem. 

15 An important step in solving i\r-dimensional assignment 

problems for N > 2 ±b finding the optimal solution of the two- 
dimensional assignment problem. In particular, we wish to 
solve the two-dimensional problem which includes the zero- 
index. Typically this problem can be thought of as trying to 

2 0 find either a minimum cost or maximum utility in assigning 
person to objects, tenants to apartments, or teachers to 
classrooms. We will follow the work of Bertsekas and call the 
first index set persons and the second index set objects. The 
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statement of this problem is given below (11.1) when the 
number of persons is m and the number of objects is n. 

Minimize X) 2^ ^a^n 

i=0 j = 0 ^ 

Subject to ^ ^17^1 i = l,...,mf 

j = 0 



_ 1 jfor all 7 = 1/ . . . / n. 

There are a couple of assumptions which we will make 
5 about (11.1), First of all, we assume that c^q and Cq-/ are 

. ? i 

7^ well-defined and the corresponding variables x^q and XQj are 

'1=^ free to be assigned- Second if a cost Cj^j happens to be 

4=1 undefined, then the corresponding variable Xj^j will be set to 

0. In effect if c^^ is undefined, we simply remove this cost 
10 and variable from inclusion in the problem. 

Notice that there are no constraints on the number of 
times person 0 and object 0 can be assigned. But notice that 
the first constraint requires that each non-zero person i must 
be assigned exactly once. Similarly, the second constraint 
15 forces each non-zero object j to be assigned exactly one time. 
Finally, the last constraint gives an integer solution, 
although we will see shortly that this constraint can be 
relaxed to admit any solution x^j-^O. One reason for requiring 
that all of the costs of the form c^q and CQj be defined is so 
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that we are guaranteed a feasible solution exists for the 
given problem. 

XI . 1 . Relaxation of One Constraint 
Relaxing the last constraint, we obtain the following 
problem: 



Minimize 



i=0 j=0 ^ j = 0 i=0 ^ 



i=0 j=0 

Subject to £ ^i-/"^! '^-^-^ i = l,.,.,zn 

j = 0 

X.^.6{0,1}. 



This lower bound is achieved and the second constraint in 
the original problem is satisfied if the following conditions 
are met : 



For i=0, 



For i^O, 



1, if c^^^u^ ^0 
0, otherwise 



1, if j = arg min{c^j^ + u^\k = 0 ^ 1 ^ ... ^ n } 
0 otherwise 



All of j^s are assigned only one time. 



The relaxation of the first constraint is analogous and would 
lead to similar results with i and j exchanged and the 
introduction of the multiplier . 

XI, 2. Relaxation of Both Constraints 
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This time we will relax both constraints and use duality 
theory to obtain a relationship between the multipliers and 
u\ 

Definition: An assignment S and a multiplier set (u^,u^) are 
said to satisfy complementary slackness - CS) if 

c.. + ul-^Uj=0 for all [i,j)ES, 
c^^ + u^^ -e for all {i,j)EA. 

Forward Auction Algorithm 

(1) Select any unassigned person i 

(2) Determine the following quantities: 

j. = arg min | c^^+ u^^| A: 6 

w.= min ^c.j^-^u^\k EAU) , k^j.^ 

In the selection of above, if a tie occurs between 0 
and any non-zero index k, then select as k. Otherwise, 
if there is a tie between two or more non-zero indices, 
the choice of is arbitrary. 

Also if A{i) consists of only one element, then set w^=^ . 

(3) Update the multipliers and the assignment: 
If ji = 0, then 
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(a) Add (i,0) to the assignment. 

(b) Update ul:--c^^. 
If ji ^ 0, then 

(a) Add Hiji) to the assignment. 

(b) Remove (i^Ji) from the assignment if was 
previously assigned. 

(c) Update u^,:= {w . ~ v ) + E = w - c., +E. 

(d) Update uj^ := '{^^ij.'^^j^ ^ -w^-E. 



Reverse Auction Algorithm. 



2 

(1) Select any unassigned object j, such that c^. + u^- >0, 



(2) Determine the following quantities: 
i^^arg min | Cj^+ u^^| A: E S ( j ) J, 

Y.= min ^Cj^,-^u^\k EB(j) , k^i.^. 

In the selection of ij above, if a tie occurs between 0 
and any non-zero index k, then select as k. Otherwise, 
if there is a tie between two or more non-zero indices, 
the choice of is arbitrary. 
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Also if B{j) consists of only one element, then set 
Yj = 

(3) Update the multipliers and the assignment: 
If ij- = 0, then 

(a) Add (0,j) to the assignment. 

(b) Update u| : = -c^.. 
If # 0, then 

(a) Add iijfj) to the assignment. 

(b) Remove [ipj") from the assignment if was 
previously assigned. 

(c) Update u-^: = u^^^ (y.-^,) +e =y.- c.^,^E . 
id) Update uj:= - ( c . + .) = "Y^.-^. 

Combined Forward/Reverse Auction Algorithm 

1. Assume that is given as an arbitrary multiplier. 

2. Adjust the value of for each object j as follows: 

If + <0, then set : = "C^oj • 

3. Run iterations of the Forward Auction Algorithm 
until all persons become assigned. 

4. Run iterations of the Reverse Auction Algorithm 
until all of the objects become assigned. 

Note at the completion of the Forward auction step we 
have the following conditions satisfied: 
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y 



11 



2 

Cq. + Uj ^0 for all objects j. 



• c.^ + u];'^Uj =0 for all S, 

• + ^ rnin|c.^ + u^^J + e for all (i,j)e S. 
Thus we can prove the following proposition. 

5 Proposition : If we assume that c^^ + ^ 0 at the start of the 
Forward Auction Algorithm and all of the persons are assigned 
via a forward step, then we have: 

c^^. + Ui+Uj ^-6 for all A. 

c^^ + Ui+u| =0 for all S. 

10 ^ij'^^j ^ for all {i,j)E S. 



Optimality of the Algorithm 

Theorem: €-CS preserved during every forward and reverse 
iteration. 

15 Theorem: If a feasible solution exists, then the resulting 
solution is with m€ of being optimal for the 
Combined Forward Reverse Algorithm. 



XII. SOME CONCLUDING COMMENTS 
2 0 Although the algorithm appears to be serial in nature, 

its primary computational requirements are almost entirely 
parallelizable . Thus parallelization is planned. 
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Step 2 is the computational part of the algorithm. 
Evaluating ^^j^-^^i and computing a subgradient use in the search 
procedure requires 99% of the computing time in the algorithm. 
This part uses two-dimensional assignment algorithms, a search 
5 over a large number of indices, and a nonsmooth optimization 
algorithm. It is the second part (the search) that consumes 
99% of the computational time and this is almost entirely 
parallelizable . Indeed, there are two-dimensional assignment 
solvers that are highly parallelizable. Thus, we need but 

10 parallelize the nonsmooth optimization solver to have a 
reasonably complete parallelization. 

If a sensitivity analysis is desired or if one is 
interested in computing several near-optimal solution's, a 
parallel processor with a few powerful processors and good 

15 communication such as on the Intel Paragon would be most 
beneficial . 

The foregoing discussion of the invention has been 
presented for purposes of illustration and description. 
Further, the description is not intended to limit the 
20 invention to the form disclosed herein. Consequently, 
variation and modification commensurate with the above 
teachings, within the skill and knowledge of the relevant art, 
are within the scope of the present invention. The embodiment 
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described hereinabove is further intended to explain the best 
mode presently known of practicing the invention and to enable 
others skilled in the art to utilize the invention as such, or 
in other embodiments, and with the various modifications 
5 required by their particular application or use of the 
invention. It is intended that the appended claims be 
construed to include alternative embodiments of the invention 
to the extent permitted by the prior art . 
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